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Cover Picture 


Two tropical weather phenomena of great Interest to meteorologists 
are prominent in this picture. Hurricane Anna (17 July 1969) revolves 
in growing intensity at the right, but the two cloud clusters in the 
center of the picture are the subjects of particular interest In Dr. 
Sikdar's paper included in this report. Perhaps cloud clusters are less 
spectacular than hurricanes, but in the long run it is probably more im- 
portant to understand their dynamics if we are to understand the meteor- 
ology of the tropics better than we do now. Most of the Sun's energy 
received by the Earth is absorbed in the tropics and cloud clusters ap- 
pear to be a key factor in the major atmospheric process of energy re- 
distribution. Since our temperate zone weather is a by-product of the 
energy redistribution process, the study of cloud clusters is a matter 
of importance. 
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Preface 


Quantitative exploitation of meteorological data from geoaynchronous 
satellites is starting to move from the laboratory to operational prac- 
tice. The papers in this report investigate several aspects of the data 
applications portion of the total meteorological satellite system. 

The first tvo papers are concerned with measurements taken from geo- 
synchronous satellite data as indicators of atmospheric physical pro- 
cesses. Professor Stremler reports in the third paper on the results of 
a preliminary study of a method to detect and measure rainfall from geo- 
synchronous orbit. While the concept does not appear feasible at this 
time, the study did Indicate that relatively modest technological ad- 
vances could provide the necessary capability in the future. The fourth 
paper documents the first half of a two-year effort to construct a 
broadly applicable model of geosynchronous orbit satellite attitude and 
pointing control systems. When completed, it is hoped that this model 
will assist systems designers and data users in minimizing and correct- 
ing geometric inaccuracies in satellite data. 

I Join with the authors in expressing sincere appreciation to the 
members of the Space Science and Engineering Center who have assisted in 
the work reported here, and in thanking the members of the National Aero- 
nautics and Space Administration who have provided the necessary sponsor- 
ship of our efforts. 


Verner E. Suomi 
Principal Investigator 
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TROPOSPHERIC VaND SHEAR AND THE RELATED SEVERE STORM CIRCULATIONS— 


A VIEW FROM A GEOSTATIONARY ALTITUDE 


B. Auvlne and D. N. Slkdar 


ABSTRACT: 

Three severe storm complexes are Investigated by means of sat- 
ellite data to reveal salient differences In circulation features, 
namely volume flux and anvil divergence. Two of these complexes 
grew In strongly sheared environments, the third In a nearly shear- 
less one. 

The analyses Indicate hardly any significant differences In 
the storm severity, growth features or circulation characteristics, 
except that In the strong sheared environment multiple cells In a 
squall-line configuration were embedded within a cloud complex 
while In the weak-shear environment only single "supercell" storms 
seemed to be favored. A possible reason for this difference may 
be due to the downward flux of momentum associated with the pres- 
ence of a high level Jet. 


1. Introduction 


Meteorological literature has concentrated most of Its attention on 
shear-related storms. Numerous case studies (e.g. Miller, 1959; Browning 
and Donaldson, 1963) have pointed to shear as an Important factor In se- 
vere storm formation. Indeed, many severe storm analysts regard vertical 
wind shear as a virtual prerequisite to severe storm formation (Fawbush 
et al., 1951; Endllch and Mancuso, 1968). 

Numerical models of convection have had varying results when shear 
Is Included In the system. Asal (1968) has Indicated that shear will sup- 
press shallow convection. Others have suggested that shear under the 
proper conditions can lead to the formation of significant and even self- 
perpetrating convection (Jlh Ping and Ll-Shoo, 1964; Takeda, 1971). 

Takeda In particular suggests that whereas steady-state storms may be 
achieved with shear as a mechanism Involved In creating or releasing 
Instability In a synoptic situation, there Is also a class of storms 
which can grow, provided the ambient shear Is very weak and the atmos- 
phere Is sufficiently statically unstable. It Is In such situations, we 
believe, that the moisture front may play a role. Carlson and Ludlum 
(1968) describe the moisture front as a mechanism capable of releasing 
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Instability In a violent way and have documented several cases In which 
synoptic conditions were similar to the case presented In this paper. 

The purpose of this paper Is to compare some of the ATS-III observed 
severe storm circulation features such as volume flux, anvil divergence 
and precipitation efficiency In a weak-shear environment (Lubbock storm 
of May 11, 1970) to the same features found in environments with strong 
upper tropospheric wind shear (In particular the cases of April 19 and 23, 
1968, and already discussed by Slkdar et al., 1970). Also briefly dis- 
cussed Is the structure of the Lubbock storm using time lapse radar data 
(WSR 57) taken from Amarillo, Texas. 


2. Synoptic Features Related to Three Severe Storms Situations 

On May 11, 1970, Lubbock, Texas was struck by two tornadoes, the 
second of which was violently destructive. The synoptic conditions that 
day Indicated a westward moving moisture front (Figure 1) along which a 
number of convective cloud masses had their origin. This moisture front 
consisted essentially of moist Gulf air flowing In to meet dry desert air, 
with little contrast In temperature between the two air masses. VAille It 
has already been documented that such moisture discontinuities In the Gulf 
region can be associated with severe convective activity (Bradbury, 1969) , 
the May 11 storms are notable for being spawned In a subtropical environ- 
ment essentially lacking In large vertical wind shear (although Reuss, 
1961, reports the existence of a giant cumulonimbus cell growing In a 
nearly calm atmosphere) . 

The situation In regard to low and high level winds on May 12 (00 Z) 
can be seen In Figures 2a, 2b while Figure 3a presents an east-west cross- 
section of Isentropes, mixing ratios and Isotachs through the convection 
area. Similar conditions prevailed 12 hours prior to the storms. At 850 mb 
one can easily see the flow of air northward from the Gulf. It Is notable, 
however, that In the moist airflow from Amarillo (AMA) to Topeka (TOP) , a 
region of relatively no shear exists as can be seen by comparing the low 
and high level wind speeds. At both AMA and TOP the wind vectors remain 
relatively constant from 30 kts at the surface to similar wind speeds at 
higher levels. Further to the south at Midland, Texas (MAP), a significant 
shear Is, however, present, due evidently to the Influence of the subtropical 
Jet over north Mexico. What Is significant In this May 11 case Is the fact 
that severe weather Is not associated with this strongly sheared regime 
where vertical differential advectlon could be expected to encourage 
convection, but with the moderate- to low-shear area further north. At 
the time of these synoptic charts, vigorous convection was to be found at 
Lubbock (about midway between MAF and AMA) northeastwards to the east of 
AMA and again In northeast Kansas. Thus, while the shear at Midland of 
about 1.7 X 10"3 sec"l Is somewhat below that prescribed by Marwltz (1972) 
for supercell storms (2.5 to 4.0 x 10~^ sec~l-), values at Lubbock likely 
averaged 0.9 x 10~^ sec~^ with even smaller values to the north of Lubbock 
and In Kansas. 

One may also note In the 200 mb map that the region of severe 
weather may be said to be In a "dlffluence" region between the polar jet 
over Nevada and the subtropical Jet over northern Mexico. A similar sit- 
uation Is to be found In other severe storm situations, as, for Instance, 
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Figure 2. b) 200 mb winds, OOOOZ May 12, 1970 with Isotachs. 






200 mb 
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Figure 3. b) croee-eection through severe storm area, 1200Z, 

April I9g 1966s 
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Figure 3. c) Same, except for OOOOZ, April 20, 1968, 








April 23, 1968. The sudden decrease In wind speeds west of northern 
Texas would seem to require a subsidence region In this same area and. 
Indeed, the area to the west of the moisture front Is perfectly clear on 
the ATS cloud photographs. How this subsidence Is In fact related. If at 
all, to the formation of convection Is not known. 

Included for comparison are the cross sections of April 19 and 23, 
1968 (Figure 3b to 3e), two days when strong shear was prominent in Its 
association with severe weather. Note that for these two days momentum 
appears to be transported downwards In the area of convection (the area 
Identified by the vertical transport of moist air above the Inversion, 
and the surface weather reports). As opposed to the May 12 storm where 
only two Isolated tornadlc storms were reported, April 19 and 23 produced 
severe weather over an area of several states. The Increase of low level 
wind speeds has been related among other factors to severe weather devel- 
opment In midlatitude frontal zones by Johnson and Sechrlst (1970) and 
Nlnomlya (1971). The May 11 case serves as a reminder by Its weak shear 
that, despite recent emphasis on the role of the Jet stream, strong shear 
In the region of convection Is not a necessary condition for severe 
weather. 

The cross-sections also reveal three different surface mechanisms 
associated with the storm outbreaks. While April 19 has a dry front, and 
April 23 a cold front. May 11 features the previously mentioned moisture 
front. The contrast of the two air masses associated with this front 
leads to conditions of high Instability with dry air having a steep lapse 
rate overlying moist Gulf air. However, an Inversion between the air 
masses discourages convection. The moist air, moving westward, has 
covered areas to Amarillo, Texas (AMA) with mixing ratios up to 12 g/kg. 
At Amarillo and roughly In a north-south line, the storms broke out late 
In the day of May 11. The Lubbock storms themselves had developed and 
were near maturity by OOZ, May 12, and are responsible for the upward 
transport of moisture at Amarillo seen on the OOZ cross-section. We may 
conclude, then, that while upper tropospheric wind shear may enhance 
growth of severe storms, some such storms nevertheless can develop with- 
out such winds, given the right conditions of instability. 


3. Data Analysis Techniques 

The first Lubbock tornado struck an hour and a half after darkness 
had fallen. The cloud photographs of storm growth from ATS-III satellite 
are available, until 0100 QIT on May 12, both for the Lubbock tornadlc 
storm and the several smaller storms around It. The elements In each of 
these photographs — five In number spaced 20 minutes apart — have been re- 
corded In terms of brightness levels on digital tapes. Mesoscale convec- 
tion systems which are difficult to analyze on photographic Images be- 
cause of their limited dynamic range can easily be treated on digital 
contour displays. Such a display gives a detailed and enlarged portrait 
of each cloud, particularly of cumulonlmbl and of their associated anvils, 
provided It Is properly brightness normalized for sun-satelllte-cloud 
geometry. The brightness of a cloud mass on a satellite photograph de- 
pends on many variables Including solar zenith angle, satellite zenith or 
viewing angle, the relative azimuth between the sun and satellite measured 
from the view spot, satellite characteristics, and the picture processing 
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Figure 4. a) ATS-III brightness contour display and photograph of storm 
clouds on 2302 Z May 11, 1970. Clouds outlined In black 
are the same as those viewed In the lower left hand corner 
of the photograph. Clouds are numbered for reference In 
flux calculations. 
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Figure 4. c) Same, except for 2346 Z. 
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Figure 4. d) Same, except for 0008 Z. April 12. 1970. 



technique. Unfortunately, to date we have no unique brightness normal- 
ization technique available. As pointed out by Martin and Suomi (1971), 
the cloud brightness largely depends on solar zenith angle, and a uniform 
thick cirrus canopy atop a deep convection system can be treated as a 
quasl-Lambertlan reflector. Although the solar zenith angle In the case 
of the Lubbock storm of May 11, 1971, Is large, we have assumed the 
visually evident uniform cirrus canopy Is such a reflector. The bright- 
ness data analyzed In this paper have therefore been corrected for sun 
angle change only, and It Is hypothesized that the time-change In the 
cumulonimbus brightness field can be attributed for the most part to 
cloud evolution process, and that the brightest portion of the growing 
cloud would correspond to active portions of the complex that had pene- 
trated to a significant depth of the troposphere and to the areas of 
precipitation (Slkdar, 1972). Four such brightness contour displays along 
with corresponding satellite photographs are presented In Figure 4. The 
clouds are arbitrarily numbered from 5 to 12. The dotted areas In the con- 
tour displays are the transposed radar echoes observed from Amarillo, Tex- 
as. In most cases a given cloud seems to be associated with only one pre- 
cipitation echo. Cloud #7 Is the tornado-producing cloud. A particular 
feature of this cloud not evident In similar studies of the April 19 and 
23 storms Is a clearly defined tower penetrating the tropopause on the 
southwest side of the cloud as evidenced by the protruded bright spot 
with Its shadow stretching to the north. A simple geometrical computa- 
tion using the length of this shadow and the zenith angle of the sun 
yields a height of about 6 km above the cirrus layer, an observation 
In close agreement with that of Roach (1967). 

As for Che determination of volume flux atop a severe storm, we have 
assumed following Newton (1966) that "essentially all the air originating 
In the updraft remains In the upper troposphere, spreading out mainly In 
the anvil plume downshear but appreciably also In the upshear side." In 
a great plain thunderstorm the cirrus shield appears several kilometers 
thick; however, only the top portion, say one kilometer, probably con- 
tains the air originating In the subcloud layer and In the remaining part 
the anvil air probably originates essentially from the strongly entrained 
air around the updraft. Fujlta and Byers (1962) show an anvil thickness 
of 1.5 kilometers from photogrammetrlc measurements of a severe storm, and 
Ludlum (1966) quotes 1.0 kilometer anvil thickness for moderate cumulo- 
nimbus clouds. In view of these observations and the fact Chat not all 
of the anvil Is part of the outflow layer, an assumption of 1.0 kilometer 
thickness for the layer of outflow seems more appropriate. In any case, 
our desire here Is to provide an approximate estimate of the actual 
volume outflow In Individual storm complexes. 

Assuming that evaporation at the anvil edge Is small (Darkow, 1963) 
and relatively constant for a 15-islnute Interval period, one may use time 
changes In anvil expansion as a legitimate tool for measuring the volume 
flux In a storm or group of storms (Slkdar et al., 1970). The anvil area 
expansion can be easily determined following normalized Isopleths of 
brightness field appropriate for cirrus shield on time-lapse digital 
displays, and then volume fluxes obtained with 1.0 kilometer thickness 
outflow layer. 

The anvil divergence can be determined from Che expression 
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A ‘ dt 


( 1 ) 


where A Is the cirrus shield area and t Is time. The magnitude of 
divergence gives a measure of the storm's severity. The volume flux can 
be related to anvil divergence by the expression 


where AZ 


dt 


A* 

H 


1.0 kilometer (assumed thickness of cirrus shield). 


( 2 ) 


4. Results: The Circulation Features as Observed from ATS-III 


Shown In Figure 5 Is a plot of volume flux versus time for the Hay 11 
storm. Numbers on the curves correspond to the numbered clouds In Figure 
4. Note that the volume flux Is not constant for all the clouds but shows 
a tendency to fluctuate. One cannot, therefore, characterize these storms 
as steady-state, at least for variations on the scale of 20 minutes. 
Similar fluctuations can be observed In the volume fluxes for April 19 
(Slkdar et al., 1970). As one might expect, the tornado storm #7 has the 
largest value of volimie flux. Although the latter observations could not 
be presented for lack of good quality digital tapes, the measurements 
obtained from enlarged photographs Indicate that the volume flux of cloud 
#7 Increased to a maximum of 1.5 km’/sec at around 0018 GMT, and there- 
after remained more or less steady until the last picture taken at 0051 
GMT. The next largest storm In terms of volume flux, cloud #8, was also 
associated with considerable radar echo activity — more so than any of the 
other neighboring clouds besides #7. 

Figure 6 gives a plot of anvil divergence as a function of time. 

While the anvil divergence magnitude atop a severe storm Is a direct 
measure of the storm Intensity (Slkdar et al., 1970), It Is not correct 
to assume that all parts of a storm cloud are equally divergent. As will 
be seen from Eq. (1), the computed divergence will depend on the size of 
the area Included In the measurement, l.e. on the scale chosen, even If 
the rate of area change may remain constant. The decline of anvil di- 
vergence magnitudes with time shown In Figure 6 Is consistent with the 
findings of Auvlne and Anderson (1972) , namely that we are dealing with 
a highly divergent core region surrounded by a region of nondlvergent 
flow resembling the behavior of a point source. The anvil divergence of 
cloud #7 In Figure 6 Is an order of magnitude higher than that of an 
ordinary thunderstorm. In agreement with Slkdar et al. (1970). One may 
conclude from these measurements that a severe storm growth need not al- 
ways be coupled directly to the presence of upper tropospheric wind shear. 

Table 1 provides a comparison of the computed volume flux and anvil 
divergence of storm #7 with those observed on April 19 and 23 In sheared 
environment. The average volume flux on May 11, 1970, was somewhat 
smaller In magnitude although the anvil divergence was significantly 
higher. Referring back to Eq. (2), one can visualize that the volume 
fluxes determined from area measurements depend on scale size of the 
clouds. The average cloud size of cloud #7 as shown In Table 1 was at 
least an order of magnitude smaller than those of April's two cases. 
Furthermore, as evidenced from the radar echo distributions, the May 11 
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storm #7 had only one solid echo while those on April 19 and 23 had 
several In them Indicating the presence of multiple point sources con- 
tributing to the average volume fluxes. Nevertheless, the volume fluxes 
and divergences of these three storm systems appear to have the same order 
of magnitude. 

On the basis of these growth parameters, one cannot distinguish be- 
tween shear-associated and shear less storms; a shearless environment Is 
capable of producing storms as Intense as any found In the midlatitude 
frontal zones. 

Severe storms have generally been found to have higher precipitation 
efficiencies than smaller and more ordinary storms. To see If the Lub- 
bock storm conformed to this expectation, we computed the precipitation 
efficiency using the mass flux at 0048Z on May 12 obtained from the last 
two photographs. The volume flux at this time had become fairly constant. 
It was assumed that the precipitation output In this storm also remained 
approximately constant during the time Lubbock, Texas, recorded precipi- 
tation (00S3Z to 041SZ) . Using this precipitation record, an Inflow 
layer Into the storm based on an LCL of 11,000 ft, an average mixing ratio 
for this layer of 9.25 g/kg, and a standard atmospheric density of .5 x 
10~3 gcm~3 at the top of the storm, one obtains a precipitation effi- 
ciency of SIX. This compares to values of lOZ for light thunderstorm 
(Braham, 1952), to values of 50X (Newton, 1966; Auer and Maurltz, 1968) 
and 60X (Fankhauser, 1971) for severe storms. 


Table 1 



Average 
cloud size 

(lOW) 

dV/dt 
(km’/ sec) 

Anvil 

divergence 

(sec“l) 

May 11, 1970 

0.15 

1.2 

1.0 X 10“^ 

April 19, 1968 

1.5 

3.2 

3.0 X lo"^ 

April 23, 1968 

2.5 

4.6 

5.0 X 10"^ 

Average for mod. 

thunderstorms 

- 

0.2 or less 

- 

5 . Some Comments 

on the Growth and 

Structure of the 

Lubbock Storm 


We have shown that from the viewpoint of the ATS satellite, the 
Lubbock thunderstorm Is, In regard to mass flux, quantitatively similar 
to shear-induced thunderstorms. We may still ask, however, whether this 
storm Is qualitatively different from Its midlatitude counterparts. In 
other words. Is the structure and growth of this storm similar to that 
of the typical Colorado and South Dakota hailstorm (documented by Dennis 
et al., 1970), or to the severe storms found In the plains states to the 
east (Browning, 1963; Newton, 1966, etc.), or Is It unlike either of 
these? 
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CYCLONIC MOTION 
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Figure 7. a) Radar (WSR-57) photographs of the Lubbock tornado echo 
(Echo I) and a nearby echo (Echo II) depicting cyclonic 
motions. Note the turning of the Echo I axis from a NE' 
SU to a N-S orientation, and the meso-low appearance of 
Echo II. 




22 



Some light can be shed on this question by looking at the Amarillo 
radar photograhs. Figure 7 reveals two Interesting features of the storms 
associated with or close to the Lubbock tornado. After 0300, Echo I, the 
Lubbock tornado cloud shows definite signs of rotary motion as the axis 
of the storm moves from a NE-SW to a N-S orientation. During the same 
period a nearby echo. Echo II, demonstrates what appears to be a meso- 
low structure such as that found by Fujlta (1958) In his examination of 
tornadlc mesoscale systems In Illinois. 

Earlier — between the occurrence of the two tornadoes — Echo I also Is 
subject to a convergence of echoes which had grown up to the rear of the 
storm. As the storm Is moving approximately northward, and the echoes 
converge In the south, growth seems to be taking place to the rear of the 
storm, the same growth mechanism described by Dennis et al. (1970) for 
South Dakota hailstorms. 

On the basis of this admittedly meager evidence, there Is reason to 
believe that the Lubbock storm resembles In growth and structure the two 
types of storms cited above, both of which are known to exist In sheared 
environments. A difference In starting mechanism and environmental wind 
conditions may not necessarily affect the structure of a severe storm 
formed. 


6. Concluding Remarks 

In this paper we have demonstrated that In regard to the nature of 
severe storm growth, namely the rapid expansion of cirrus shield In ad- 
vance of the tornado occurrencea, the Lubbock storm closely resembles the 
two severe storms of April 19 and 23, 1968 (Slkdar et al., 1970) and that 
the magnitudes of anvil divergences and mass fluxes are comparable. In 
addition, several structural characteristics cited by Investigators re- 
lating to shear affected storms seem to apply to the Lubbock storm as 
well. These characteristics Include the rotating cells and the growth 
of the cells by means of feeder clouds to the rear of the storm. Such 
observations, however, are not conclusive and further study will be 
needed to determine the exact nature of circulations In nearly shearless 
storms In contrast with what has been learned of severe storms experi- 
enced In midlatitude frontal zones. 

The difference In starting mechanism In the subcloud layer and en- 
vironmental wind conditions as discussed In Section 2 may not necessarily 
affect the growth rate of a single Isolated severe storm cell. However, 
such factors may be Importantly related to the occurrences of widespread 
multiple convection cells as seen In the storm complexes of April 19 and 
23, and their nonoccurrence on May 11. Whereas In the former cases the 
presence of a Jet stream In the region of convection allows for the 
formation of a large region of synoptic scale Instability and the subse- 
quent formation of a large, severe, multi-storm complex, the lack of such 
a condition In the Lubbock region allowed only the formation of one 
large "super-cell." 

That the presence of wind shear In some way favors the formation of 
larger or more Intense storm complexes Is Illustrated In Figure 8 where 
cirrus shield area expansion rate Is plotted against the maximum wind 
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In cloud layers (850-200 mb) In the vicinity of these three storms. Of 
course, the most Intense outbreak of severe storms Is associated with high 
wind shear while the Isolated May 11 storm has hardly any. The mechanism 
of Interaction between the severe storm growth and the environmental wind 
shear Is not yet clearly understood; It may be a result of a "chimney ef- 
fect" where the flow past the top of the rising updraft air tends to In- 
crease Che rate of mass transport through the storm. The transport of 
environmental momentum downward as evidenced In Figure 3 may also encour- 
age the formation of a large area of severe convection. Time lapse Infor- 
mation from ATS on the evolution of such severe storm complexes on a 
shorter time-scale should hopefully reflect on the multi-cell growth 
mechanism In the near future. 
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ON SOME KINEMATIC PROPERTIES OF THE TROPICAL ATMOSPHERE 


AS DERIVED FROM CLOUD MOTION VECTORS 


D. N. Slkdar 


ABSTRACT: 

Using the WINDCO system (Smith and Phillips 1972), vectors of 
small cloud movement were obtained from ATS-III picture sequences 
for the period 26-28 July 1969 In the tropical Atlantic. These 
were analyzed for various kinematic properties such as cloud tra- 
jectories, divergence, relative vortlclty and specific kinetic 
energy. Among the significant findings are the following: 

1. A significant latitudinal wind shear exists north and south 
of the ITCZ, especially In the zonal component. 

2. Cloud clusters are found In the general areas of converg- 
ence (-D) and positive relative vortlclty (C) . 

3. The vortlclty magnitudes appear to be scale-dependent while the 
divergences are not, at least In the scale 100-600 nautical miles 

on a side. 

4. The grid point values of C and D In large nephsystems 
are approximately related by the empirical equation C “ -1.5 D 
Implying that the low level divergence decreases with the Increase 
of relative vortlclty In an active cloud cluster; however, for 
some cloud clusters this relationship Is missing. 

5. Specific kinetic energy fields for the three days Investi- 
gated show regions of maxima at the locations of strong zonal 
wind shear. 

It may be concluded from this preliminary Investigation that the 
WINDCO data are very useful for examining various scales of tropical 
circulations, especially In the data-vold regions, and has good 
potential for input to a global numerical prediction model. 
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1. Introduction 


A recent technological development (WINDCO: an electronic animation 
system for obtaining accurate cloud motion) at the Space Science and En- 
gineering Center (SSEC), University of Wisconsin — Madison, has made pos- 
sible the measurement of cloud displacements to an accuracy of less than 
3 knots, a primary requirement of the Global Atmospheric Research Pro- 
gram (GARP) . Attempts are in progress to employ these data to Improve 
understanding of atmospheric circulations, especially in the data-vold 
tropical oceanic regions. 

The visible sensors of the geostationary satellites ATS-I and III 
have limited resolution in the vertical. However, from a comparison of 
conventional wind data with small cloud motion vectors, Hasler (1971) 
found the best fit level for low cloud trajectories to be around 950 mb, 
which approximately corresponds to the base of tropical cumuli. For the 
950 mb level, 68Z of the direction differences were less than 30°. 

Fujlta et al. (1969) also stated 950 mb as the best fit level while 
Hubert and Whitney (1972) and Serebreny et al. (1969) arrive at 950- 
850 mb from their comparisons. 

Future geostationary meteorological satellites equipped with Infra- 
red sensors will greatly reduce the present uncertainties of cloud height 
determination of large and thick clouds, but for moving thin cirrus 
shields the uncertainties due mainly to emlsslvlty and transmissivity 
will remain. In this connection It should be mentioned that a prelim- 
inary Investigation In progress at this Center suggests a reasonable 
correlation between the cloud brightness and thickness (Park 1972). If 
this Is so, one should be able to assign a cloud displacement to a level 
depending on the Integrated cloud brightness. Pending establishment of 
such a relationship, we will proceed with the hypothesis that small 
clouds In ATS-III Images are at the trade cumulus level. 

The need to Identify small clouds on a series of successive pictures 
at Intervals of 23 minutes Imposes some restrictions on the data selec- 
tion scheme. Nevertheless, our experiences are that the numbers of even 
such restricted cloud displacement vectors are enough to give adequate 
two-dimensional coverage of kinematic parameters such as stream function, 
vortlclty and divergences. Of course, one should keep In mind chat such 
data on occasions may be very few In areas dominated by subsidence as, 
for example. In the subtropical oceanic highs and ridges. 

The purpose of this paper Is to demonstrate Che potential of WINDCO- 
derlved cloud motion vectors for synoptic and dynamical analyses of 
tropical weather systems In the absence of upper air data. 


2. Data Selection and the Analysis Technique 

During the development and testing phase of WINDCO, cloud displace- 
ment vectors were generated for three days — 26, 27 and 28 July 1969 — 
from Phase IV of BCMEX. These data were quality controlled and access- 
ible, and therefore were selected for analysis of various kinematic 
properties. 
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Figur* la. Satellite photos shoving cloud cluster activity In the 
western Atlantic on 26 July 1969. 
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Figure lb. Same as la except for 27 July 1969. 
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Figure Ic. Sane as la except for 28 July 1969. 
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The WINDCO system measures cloud motions by tracIcioK clouds over two 
digital pictures and applying a two-dimensional cross-correlation analy- 
sis. A precise alignment of ATS picture sequences Is crucial to the ac- 
curate measurement of these cloud displacements. The second Important 
step Involves finding the coordinates of each selected cloud tracer on 
Che digital pictures, the primary criteria being that the cloud tracer 
should be found In all the pictures In the selected sequence and Its lo- 
cation should not be near the limb. For details the reader Is referred 
to Smith and Phillips (1972). 

The cloud displacement vectors produced by WINDCO were first resolved 
Into u and v components and averaged In an area of 2.5* x 2.5°. This 
procedure gave a matrix of uniform field of u and v In Che latitude 
band 10° - 30°N and longitude band 50° - 90°W. Approximately 150 vectors 
were obtained on the 26th and 28th with over 400 vectors on the 27th In 
the area bounded by 5°S - 40°N and 10° - 100°W. The time Intervals be- 
tween picture pairs on Che 26th and 28th were 26 minutes In all four cases 
but on the 27th the first Interval was 49 minutes and the second 38 
minutes. Doubtful data, namely, wind vectors with speed greater than 
30 m/s or direction opposed to smooched stream line flow (+ 20°), were 
considered Co be mlstracks of clouds or at different heights than the 
majority of Che winds and therefore were excluded from the analyses. 

The number of such rejected vectors was hardly 5-6Z. 

A contour analysis of u and v components derived from cloud dis- 
placement vectors revealed the horizontal gradient field from which basic 
meteorological parameters such as divergence and relative vortlclty fields 
were constructed as a function of time. In addition, the cloud displace- 
ment vectors were averaged to study their latitudinal and longitudinal 
variations. Also produced from the u and v fields were specific 
kinetic energy fields as a function of time. 


3. Satellite Time-Lapse Pictures and Cloud Motion Vectors 


A Satellite View of the Cloud Clusters on 26, 27 and 28 July 1969 

A cloud cluster Identified as A In Figure la at 57°W, 15°N In the 
western Atlantic was the only active cluster In the selected grid zone 
on 26 July 1969. Associated cloudiness appears to extend southwestward 
over South America. The cluster moved westward to the 63°W meridian on 
27 July, Implying a speed of ^ 5°/day and became disorganized as it 
passed to the east of Cuba on 28 July. Initial westward movement be- 
came northwestward on 28 July. The cloud cluster B located at 54°U, 

1S°N In Figure lb was just entering the grid zone In Figure la and was 
apparently linked to the dominant cyclonic circulation to Its east- 
southeast at 35°W, 10°N (tropical storm Anna). A wide field of stratocumulus 
clouds Is seen to the north of this latter system extending to latitudes 
as high as 25°N. The cloud cluster B can be followed to the location 
62°W, 15°N on 28 July (Figure Ic) thereby Indicating an apparent speed 
of 8° of longitude per day, which Is faster than cluster A on 26-27 July. 

The Implication that the trade-circulation was stronger on the 27th In 
the easter part of the zone than on the Z6th Is confirmed by the fields 
of cloud motion vectors discussed below. 
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Figure 2a. Cloud trajectories obtained from WINDCO on 26 July 1969 in 
the western Atlantic area. 




Figure 2b. Same as 2a except for 27 July 1969 and In the whole Atlantic. 





Figure 2c. Same ae 2a except for 28 July 1969. 



Circulation Features as Revealed From the Cloud Trajectories 


Figures 2a, b and c present the corrected cloud motion vectors obtain- 
ed from WINDCO. One can see from these diagrams that the data density Is 
high, much higher than once can expect from land based rawln network. 

For comparison Che streamline flow charts at 3000 ft obtained from Che 
WMO Regional Center at Miami are presented In Figures 3a, b and c. The 
subjectivity Introduced In the analysis of widely spaced rawln, aircraft 
and ship data In these charts Is evident, A comparison of Figures 2 and 
3 reveals In general an agreement between 3000 ft level streamlines and 
the low cloud trajectories presented. 

Cloud trajectories on the 26th Indicate a cyclonic circulation 
centered near 47°W, 11°N. Two dominant subtropical antlcyclonlc vortices 
are located at around 50°W and 80°H. Cloud motions were generally from 
east to southeast In the eastern half of the box; In the western sector 
they were light to variable. A pronounced latitudinal wind shear can be 
seen In the Caribbean Sea area to the north of the ITCZ at 12*N but to 
the west of the active cloud band. 

On 27 July (Figure 2b) , there appears a change In the circulation 
regime; northerly to northeasterly flow seems to dominate south of 12*N. 

A quasi-steady northeasterly flow Is found In the eastern half of the 
extended grid box while winds appear light and variable In the western 
Atlantic. Two cyclonic vortices, one at 15°N, 25*W and the other at 10°N, 
40*W, are evident on this chart. The antlcyclonlc vortex observed north 
of Cuba on the 26th moved a few degrees westward, A weak trough aligned 
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Figure 3a. Conventional wind analysis at 3000 ft level In the area shown 
In Figure 2a for 26 July 1969. 
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rigur* 3b. Same ai 3a except for 27 July 1969 and over the whole Atlantic. 
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Figure 3c. Same as 3a except for 28 July 1969. 


northwest-southeast Is located In the general area extending 15°N through 
5°N along the 67°W meridian. An eastward moving midlatitude cold front 
extending as low as 27*N In the western Atlantic was captured by expand- 
ing the data coverage northward. With only data from conventional sys- 
tems over the ocean, It Is not uncommon to miss such a midlatitude trough. 
Strong latitudinal wind shear persists In the Caribbean Sea east of 
Panama. 

By 28 July the wave to the north of the ITCZ moved to the 72“W meri- 
dian. The antlcyclonlc vortices to the northeast and northwest also 
moved westward a few degrees. A new cyclonic circulation Is seen 
centered approximately around 7°N, 78*W. 

The details In these wind fields offer a possibility of testing the 
"easterly wave" concept (Palmer, 1952; Rlehl, 1954), which holds that 
wave perturbations of the easterly trade flow are associated with 
disturbed, showery weather of the kind to be expected under a cloud clus- 
ter and that this disturbed weather typically occurs on the upwind side 
of the wave trough. To examine whether the cloud clusters are related 
to any such wind perturbations, the low cloud trajectories were carefully 
scanned and approximate locations of the clusters circled. Offhand, one 
can hardly Identify any clear-cut wind discontinuity In the cloud cluster 
area except on a few occasions; for example, at the location of cluster 
"a" In Figure 2b where a very weak amplitude wave Is Indicated. In con- 
trast to this, one would find an area of weak convergence near cloud 
cluster "E" marked by relatively stronger wind In the upwind than In the 
downwind from the cluster. In general, however, one finds a stronger 
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Figure 4b. Seme as 4a except for 27 July 1969. 
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trade flow north of a cloud cluster location than to the south and that 
all clusters may not be associated with wave troughs. 

These cloud trajectory charts clearly demonstrate superiority In data 
density over the oceanic regime where. In spite of the uncertainties In- 
volved In the estimation of cloud heights, an analyst Is able to follow 
the cloud- cluster-scale circulation more closely than otherwise would be 
possible. In most areas the data density satisfies the two-degree hori- 
zontal grid mesh spacing considered necessary for a global numerical 
model. 


4. Mean Zonal and Meridional Winds 


From the Isotach analyses of zonal and meridional components of cloud 
motion vectors, means as a function of time, latitude and longitude are 
constructed (see Figures 4a,b,c). The number of observations used In 
finding the means are Indicated against each data point. A significant 
zonal wind shear north and south of the ITCZ at 12*N Is apparent In these 
diagrams especially on the 26th and the 27th. Although not very pro- 
nounced, the meridional components exhibit confluence to the south of the 
ITCZ, and particularly so on the 27th and 28th of July. 

The longitudinal distributions of mean u and v fields show pat- 
terns seemingly Influenced by the location of vortices. In the zonal 
profiles the low magnitude between 70° - 75°W and also around 85 - 90°W 
meridians may be related to a col region between the two antlcyclonlc 
vortices. Interestingly the locations of these regions of low zonal 
wind speed are constant on the three days Investigated In spite of slow 
westward movement of the antlcyclonlc cells. While the zonal winds ex- 
hibited more or less the same pattern, the meridional wind distribution 
significantly changed; namely, the bump Indicating southerly flow of 
3-4 mps at around 63°W on 26 July (Figure 3a) became weaker on subse- 
quent days and shifted westward to approximately 72°W on 28 July. In 
the vicinity of cloud clusters the average zonal shear was found 
= 7 X 10"6 aec”^ (cyclonic) and the average meridional shear was around 
3 X sec~l (cyclonic). 


5. Divergence and Relative Vortlclty Fields Related to Cloud Population 


Divergence Fields 

Figures 5a,b,c are low level divergence fields produced from cloud 
motion fields upon which are superimposed clouds from ATS photos. The 
major areas of convergence are found located near the cloud cluster In 
agreement with what Hasler (1971) has found. The magnitudes of con- 
vergence ranged from 3 x 10“5 sec“^ to 9.0 x 10~5 sec“^ depending on 
the Intensity of cluster activity. In Intense growing storms the magni- 
tudes may be as high as 10 ~^ sec~l (In 2 ° x 2 * zone) which Is about two orders 
of magnitude higher than the value 10"4 sec"! frequently quoted of large scale 
motions. The divergence magnitudes In the clear environment ranged 
between 0.2 - 5.0 x 10~^ sec"! on the average. There are occasions 
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Figure 5c. Same as 5a except for 28 July 1969. 
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when clouds are present In the areas of no significant convergence field. 
This may be attributed to Inactivity within the cloud clusters at the 
time of observation; some might simply have been composed of debris left 
over from previous convection. To test this suggestion one could check 
the cirrus expansion rate atop a cluster using Slkdar's (1969) technique. 

While editing and smoothing Is necessary because some signal noise 
Is unavoidable, the extent to which this can be done without concealing 
desired meteorological Information Is yet an unanswered question. Space- 
smoothing Is effectively done by the choice of grid size. It Is generally 
accepted that meteorological parameters like divergence and vortlclty 
fields are scale dependent. To test the validity of this concept we 
have plotted In Figure 6 divergences obtained from movleloop technique 
after Hasler (1971) averaged In 5* x 5‘ and 10' x 10' (latitude x longi- 
tude) grid spacing for the same area and date against those obtained 
from the WINDCO system (spacing 2.5' x2.5' on the 26th and 28th while 
2' X 2' on the 27th). The absence of any large scatter probably suggests 
that the low level divergence or convergence around a tropical cluster 
hardly depends on the scale size In the range 100 to 600 nautical miles 
on a side; In other words, B and C scales are effectively coupled to 
large scale motions. In Interpreting these data one should realize, 
however, that some discrepancies will be obvious because of the differ- 
ences In the technique and also because the same cloud may not have been 
followed as a tracer by the two techniques. Even so, the Independently 
derived wind fields exhibited excellent agreement (Hasler, 1971). 


Relative Vortlclty Fields 

Figures 7a,b,c show clouds superimposed on relative vortlclty fields 
for the three days Investigated. Patterns are complex and occasionally 
difficult to Interpret. In general the large clouds are located westward 
or downstream In the regions of large positive relative vortlclty with a 
vortlclty maximum often located to the west of the cloud cluster. The 
large positive values to the east and northeast of the cloud cluster are 
probably related to wind shear. 

Figure 8 Is a scatter diagram Illustrating the effect of scale size 
on the averaged relative vortlclty fields. A wide scatter probably 
means that the relative vortlclty field Is more difficult to predict from 
large-scale motions than the divergence field and Is scale dependent. A 
spatial analysis of components 3u/3y and 3u/3x (not presented here) 
reveals that a major contribution to relative vortlclty field In the 
cluster area comes from the shear In the zonal flow, l.e., - 3u/3y . 

The ratio of the zonal to the meridional shear ranges from 2:1 to 5:1 de- 
pending on the Intensity of cloud activity. While In the vicinity of a 
cloud cluster the shear Is cyclonic; the clear environment exhibits 
predoDolnantly antlcyclonlc wind shear. 


Relationship Between Divergence and Vortlclty Fields 

In this section we look for relationships between the divergence and 
vortlclty fields related to the location of a cloud cluster. Of the 
three days presented, the 27th has the best data coverage and two active 
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clusters In the central North Atlantic, with an Intense cyclonic circula- 
tion near the ITCZ In the eastern North Atlantic. Therefore, only the 
analysis for the 27th Is presented. 

Figure 9 Indicates the locations of maxima of relative vortlcy (C) 
and of divergence (D) related to the cloud cluster. One clearly finds 
In this diagram that the locations of and are separated 

by a considerable distance. In the case of Intense cyclonic circula- 
tion (Anna), however, both (cyclonic) and (convergence) are 

located near the storm's circulation center; this may Imply a different 
flow pattern than Is typically associated with a cloud cluster. 

Although the locations of C and -D and the cloud cluster 
° max max 

do not exactly coincide, the cloud cluster growth Is closely dependent 
on the convergence, for the latter becomes a dynamic necessity for the 
air to rise to the cirrus outflow layers through convective towers 
embedded In active cloud clusters. Charney and Ellassen (1964), 

Ooyama (1964) had proposed CISK (Conditional Instability of the Second 
Kind) as an Important mechanism for the Intensification of cloud clusters 
and tropical storms. If CISK Is the mechanism of convergence, one would 
expect a clear-cut relationship between the grid point values of C and 
D obtained from this analysis, and the Intensity of cluster activity. 

We assumed here that the small clouds used as tracers of air motion are 
tied to the boundary layer dynamics. 

Figures 10a and b are scatter diagrams, presenting grid point values 
of C and D Inside and near the cloud edges. The diagram reveals In 
general that In the cloud cluster vicinity, the divergence or convergence 
usually do not exceed 2.5 x 10“5 sec“l and the relative vortlclty Is pre- 
dominantly antlcyclonlc; the cloud cluster area la characterized by 
positive relative vortlclty and negative divergence and Che grid-point 
values of C and D Inside an Intense cyclonic circulation (for example, 
cloud #1 In Figure 9) are located In the same domain. The magnitudes 
are, however, significantly higher than In an ordinary cluster, and the 
relative vortlclty Is dominant. An approximate relationship 
C ” -1.5D emerges from the scatter diagram Implying a decrease of low 
level divergence with Increasing relative vortlclty In the ratio l:-0.67. 


7. Specific Kinetic Energy 


Figures 11a, b, c present specific kinetic energy fields defined by 
7(u^ -t- v^) for 26, 27 and 28 July, respectively. No specific conclusions 
can be drawn from such an analysis except that one finds consistency In 
the location of specific kinetic energy generation sources on the three 
successive days Investigated. The maxima are approximately where horl- 
zonal wind shears have the highest values, l.e., at the southern 
periphery of the antlcyclonlc vortex In the western as well as In the 
eastern North Atlantic (Figures 2a,b,c). Two dominant source regions are 
noted during this study period: one centered a:round 17*N, 60*W and the 

other around 15*N, 80°U. On extended charts of the 27th, another source 
region Is around 12*N, 40*W on the ITCZ west of Africa. 
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Figure 11a. Specific kinetic energy field constructed from 26 July 1969. 
data. 
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Figure 11b. Same as 11a except for 27 July 1969. 






Figure lie. Same as 11a except for 28 July 1969. 
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8. Concluding Remarks 

In this sample study of the one- layer model, we have seen that the 
low level cloud trajectories obtained from the VINDO system exhibit pat- 
terns comparable to streamlines at lower tropospheric levels (3000 ft). 

In the u and v fields analyzed from these data we do see signifi- 
cant changes In large-scale motion field on a time scale of a few hours. 
Assuming the wind finding technique to be working proficiently, this 
short time change may be attributed to cloud evolution processes or real 
changes In the horizontal wind field, for the cloud targets selected as 
tracers may not be Individual cumuli but small clusters and patches. 

Also seen In the u and v fields Is a significant latitudinal wind 
shear north and south of the ITCZ, especially In the zonal wind. 
Meridional components also exhibit shear south of 1S*N. 

Superimposing divergence and vortlclty fields on the cloud field, 
one gets the general Impression that the cloud clusters are In the 
areas of convergence and positive relative vortlclty. The vortlclty 
magnitudes appear to be scale-dependent while the divergences are not, 
at least In the scale 100-600 nautical miles on a side. The grid point 
values of C and D are found approximately related by the equation 
C ■ -1.5D, Implying that the low level divergence generally decreases 
with Increasing relative vortlclty In an active cloud cluster; however, 
for some clouds such a relationship Is missing. One may conclude from 
this sample study that the CISK mechanism may be valid for large tropical 
nephsystems with Intense convective activity In them — but may not be true 
for all tropical cloud clusters. 

Finally, the specific kinetic energy patterns consistently show two 
dominant source regions to the north of the ITCZ during the study period 
which are possibly related to zonal wind shear. 
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THE FEASIBILITY OF THE APPLICATION OF A GEOSTATIONARY SATELLITE 


RAKE SYSTEM TO MEASUREMENTS OF RAINFALL 


F. G. Stremler 


This Is a report on a feasibility study of the use of a RAKE radar 
system In a geostationary earth satellite to measure rainfall budgets for 
large, normally Inaccessible regions of Che earth's surface. A subcon- 
tract report of the Collins Radio Company Is appended as a major portion 
of this study [1]. Professor Blrkemeler of the Department of Electrical 
Engineering has been very successful In applying RAKE techniques to tropo- 
scatCer measurements [2], and his comnents and advice have also been 
valuable In this study. 

Major problems confronting Che use of a RAKE radar system to synchro- 
nous meteorological coverage are: 

(1) a doppler spectrum, arising from rainfall, that Is not sharply 
defined (l.e., large variance); 

(2) an unfavorable slgnal-to-clutter ratio resulting from the rela- 
tively small returned radar power from rainfall In the presence 
of a strong return from the earth; and 

(3) lack of a suitable satellite for a RAKE experiment as a result 
of the high power and frequencies. 

These are discussed here with particular reference to the Collins report. 

Because of Che distribution of drop sizes In precipitation, a radar 
pointing vertically will observe a doppler spread resulting from a dis- 
tribution of fall velocities. Measurements and data In the literature 
on doppler spectra from rainfall are very meager. Calculations made by 
Collins (see their Figure A-6) show that the standard deviation of fall 
velocities Is about 2 m/sec and Is fairly Independent of rain Intensity. 
Calculations using a Marshall-Palmer distribution for drop size [3] and 
measurements of fall velocities versus drop size found In Medhurst [A] 
yield very similar results, with slightly smaller variance. Lhermlte 
aa quoted by Nathanson [5] , states Chat Che standard deviation of fall 
velocities Is approximately 1.0 m/sec; Nathanson also adds a turbulence 
factor of about 0.7 m/sec. These references and a limited set of 
measurements reported by Nathanson and Davidson [5] reveal that the 
assumption of a standard deviation of 2 m/sec, as used by Collins, Is 
certainly not unreasonable though perhaps a little pessimistic. The 


62 


outcome of this Is that even If the earth return were not present, the 
minimum doppler bandwidth arising from rainfall Is on the order of 
400 cm/sec or ~ 2v/X - 800/X Hz where X Is In cm. An earlier 
estimate of performance by Collins [6] was based on the possibility of 
using a doppler bandwidth Bj - 1 Hz resulting In significantly better 
receiver sensitivities. 

It Is assumed that the tropics Is the region of Interest for RAKE 
rainfall measurements. In this region the source of heavy rainfall Is 
primarily from cumulonimbus clouds, and the water content and charac- 
teristics of such clouds are fairly well known [7], [8]. Representative 
numbers chosen In Che Collins report to compute signal attenuations are 
not extreme and appear to be reasonable for a 95% probablllty-of- 
detectlon performance [9]. Expected radar signal return levels calcu- 
lated on the basis of a transmitted power of 100 watts and the attenua- 
tions used by Collins are shown In Figure 1 together with receiver noise 
levels corresponding to a 1 Hz and a 800/X Hz doppler bandwidth. As Che 
transmitted power levels are changed, the returned power levels laove on 
the graph but the receiver noise thresholds remain fixed. Conclusions 
are that about ISO watts are needed to detect a rainfall rate of 10 nn/hr 
over a region 50 km In diameter at X - 2 cm using a 10 m antenna 
(Ideal) (also see Figure A-10 In Che Collins report). As a more practi- 
cal case, 1 kw Is required at a wavelength of 4 cm If the antenna gain Is 
limited to so dB. These numbers are determined by adding the required dB 
Co raise the rainfall return power curves In Figure 1 to match Che re- 
ceiver noise levels. 

To complicate matters, a very large earth return arrives at the re- 
ceiver a few microseconds after the rainfall return. This received power 
level (PrE) large compared to the earth return, as Illustrated In 
Figure 1. There are several options open to suppress this strong ground 
clutter return. The first method assumes chat the doppler of the satel- 
lite with respect to the earth can be removed and that the earth return 
can then be separated from the rainfall return on the basis of doppler 
shift. This Is Che assumption used In Che Collins report. 

Doppler clutter return for land masses will be very low while doppler 
returns from waves will vary with sea state, but will rarely exceed 
2 m/sec [10]. Therefore, Che doppler spectra from Che earth return and 
the rainfall return are nonoverlapping and can be separated by filtering. 
Assuming a 40 dB attenuation on this basis, the remainder of the attenua- 
tion required (on the order of SO dB) must be picked up using the proper 
signal design. Collins used the clutter rejection capabilities of a RAKE 
signal long enough to fill the time delay to the earth and return. 

Several alternatives can be considered In the earth return problem. 

One method Is to use short coherent code words spaced just far enough 
apart so that the earth return Is received between words. Another possi- 
bility Is to employ frequency agility Instead of time agility In separat- 
ing the rainfall returns from Che earth returns. Answers to how these 
signal changes affect RAKE system performance are not presently known [11] 
and furnish an Impetus to new research In signal design In this area. 

These methods, while perhaps promising better clutter rejection. Increase 
the complexity of the satellite RAKE system and should be studied in depth 
before dependence on them Is assumed. 
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Figure 1 
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A third area of Investigation Is the availability of existing or 
planned satellites on idilch a RAKE rainfall detection experiment could be 
attempted. The Collins report cites the status of the ATS-G. Availabil- 
ity of a suitable satellite or pair of satellites Is generally precluded 
for one or more of the following reasons [12]: 1) Radio frequencies above 

6 GHz are not currently used In satellites (except ATS); 2) few satel- 
lites, If any, could meet the power and/or bandwidth requirements needed; 
3) satellite communication systems frequencies are chosen so that no fre- 
quency used on the downlink of one satellite Is used on the uplink of 
another. Some military satellites might have the capability to try such 
an experiment but data have not been available on this. 

In conclusion, the Collins report has correctly defined the problem 
and estimated the feasibility of the RAKE radar rainfall measurement 
based on available theory and data from the literature. This study points 
out several main problems In the Implementation of such a system and areas 
for further work and Investigation: 

(1) Can the doppler dependence on rainfall rate be checked with 
experimental data? 

(2) Given these data, can rainfall rates be reliably made using 
doppler estimators? 

(3) Can suppression of earth return clutter be Improved using dif- 
ferent signal design techniques? 

(4) Can a satellite be found to try the experiment without excessive 
changes In the satellite? 
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APPENDIX 

The Feasibility of the Application of a Geostationary Satellite RAKE 
System to Measurements of Rainfall* 

1. Introduction 


Preliminary analysis has suggested the use of the unique radar mapping 
capabilities of a RAKE system Incorporated Into a geostationary satellite 
Implementation, as graphically Illustrated In Figure A-1, to obtain mea- 
sures of rainfall budget for the relatively large, but normally Inaccess- 
ible, tropic regions of the earth. 

The geostationary satellite RAKE approach to measurement of rainfall 
budget Involves the transmission of a probing signal consisting of a 
microwave frequency carrier, biphase modulated by a bilevel pseudorandom 
noise sequence, from the satellite toward the earth with vertical or near- 
vertical Incidence to the earth's surface. Portions of the transmitted 
signal are reflected from rain at various altitudes and from the earth's 
surface to return to the satellite as a signal similar to the transmitted 
signal but dispersed In both time and frequency dimensions. The satellite 
directs this dispersed signal to the RAKE receiver wherein the signal Is 
cross-correlated with appropriately delayed replicas of the transmitted 
signal to effectively separate the received signal Into a series of time 
delay Increments. Since the desired rain return signal, which corresponds 
to rain observed In a given altitude region, differs In delay from rain 


^Prepared by R. H. Pool, T. L. Rise and A. R. Hamilton; Collins Radio 
Company, Cedar Rapids, Iowa, August 1972. 
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RAKE TRANSMITTER AND RECEIVER 



Figure A-1. Geostationary satellite RAKE rainfall measurement concept. 



returns of other altitude regions and from the earth return, the desired 
rain return signal can Ideally be Isolated from the others through the 
RAKE receiver cross-correlation process. The recovered rain return sig- 
nal can then be spectrally processed to obtain a doppler frequency spec- 
trum that relatea to the distribution of rainfall velocity within the 
observed altitude region. From the relationship of rainfall velocity 
and raindrop radii dlatrlbutlons and a knowledge of the antenna beam- 
filling factor of the observed rainfall region, perhaps determined op- 
tically, the rain return signal power and Its frequency distribution can 
be employed to obtain a measure of rainfall rate. Subsequent time Inte- 
gration of rainfall rate will allow a value for total rainfall to be 
determined. 

This report relates the procedures followed and the results obtained 
during the performance of an In-depth analysis of the above described 
geostationary satellite RAKE system concept to determine the feasibility 
of Its application to measurement of rainfall budget. Consideration Is 
given to optimization of the system concept, determination of a "best" 
technique for Implementation and the possibilities of Implementing the 
system on the basis of planned geostationary satellite configurations. 


2. System Analysis and Optimization 

In this section the basic equations, parameters and parameter relation- 
ships defining the geostationary satellite RAKE system are presented. 

These are then applied to the problems of optimization of the system con- 
cept on a theoretical basis without regard to a specific form of system 
Implementation, although practical limitations and constraints are In- 
cluded In the optimization procedurea. 


2.1 Basic Equations and Parameters 

Two equations relating the signal power returned to the satellite from 
rain and the earth's surface form the basis for determination of system 
feasibility. These two equations Involve relationships of the many 
meteorological and radio parameters that define the characteristics of 
the radio propagation path and terminal equipment. 

The first of these equations, basically from Battan [1], defines the 
total average signal power returned to the satellite from spherical water 
particles contained In an observation region above the earth's surface 
having a volume determined by the area Illuminated by the antenna and a 
height h. 


P 


raR 


P^A^eWFlKin ^ 
” X* r* 



10 


-0.2L 


( 2 - 1 ) 


where 
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“ Average rain return signal power (watts) 

- Transmitted power (watts) 

• Apertural area of satellite antenna (cm^) 

6 - Antenna azimuth beamwldth (radians) 

^ - Antenna elevation beamwldth (radians) 

h •• Height of rain observation zone (cm) 
ifi m Beam-filling factor (Z/100) 

F - Correction factor (see L. J. Battan, 1959) 

|k|^ - Complex Index of refraction (see L. J. Battan, 1959) 

X - Wavelength of transmitted center frequency (cm) 

Z ” Reflectivity factor (cm*/cm’) 

r - Range (one-way, satelllte-to-ref lector distance) (cm) 

k - Propagation attenuation rate, one-way (dB/km) - k -t- k -I- k 
(gases, k^; clouds, k^; precipitation, k^) * ** 

L - Equipment losses, one-way (dB) 


The second equation, based upon work by Beckmann and Splzzlchlno [2], 
defines Che received power reflected from a diffuse scattering sea sur- 
face. 


P A*e<(ilO“'‘ 
t P 

36X*r* 


cot^Bo exp[-tan*a/tan*Bo] 


where P^.^ Is the received earth return power, a Is the angle of signal 
Incidence to the earth's surface measured from vertical, and Bo Is the 
mean slope of the sea surface Irregularities. The remaining parameters 
are as previously defined. 

These two equations will be used In various forms throughout Che re- 
mainder of the report. The need for the rain return power equation Is 
perhaps obvious. The earth return power equation enters Into the prob- 
lem due to the finite clutter rejection capability of the RAKE system. 

If the earth return power exceeds the rain return power by an excessive 
amount, the correlation processes of the RAKE receiver and the spectral 
processing Chat follows may not be capable of adequately Isolating the 
two signals. The sea surface return, much stronger than the ground 
return, provides for a worstcase, but certainly reasonable, situation. 
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2.2 Parameter Relatlonahlps 

The progression toward system optimization and finally the determin- 
ation of system feasibility on a theoretical basis requires a luiowledge 
of the parameter relationships presented In the following paragraphs. 

Since the basic power return equations expressed In the previous para- 
graph 2.1 are functions of such a large number of variables. It Is neces- 
sary to develop these parameter relationships In order to reduce the num- 
ber of variables Involved In the basic equations. Also, reasonable con- 
stant values will be chosen for certain parameters. Further, the range 
of certain variables will be restricted, from practical considerations, 
to provide additional simplification relative to the totally general case. 
After optimization Is completed. It will be possible to evaluate the ef- 
fects of variations In the restricted variables as perturbations rather 
than true variables. 

2.2.1 — Pseudorandom Sequence Clock Frequency . Starting with a given 
value for the height of the rain observation region h, the effective 
transmission bandwidth W can be determined as 



(2-3) 


where c Is the magnitude of the free-space propagation velocity. Us- 
ing the definition of W provided by Berkowitz [3], W Is found to also 
be approximately equal to the pseudorandom sequence generator clock fre- 
quency f(. If the actual RAKE transmitted signal bandwidth Is 21^, 

With fc set equal to W, then f^ Is defined as a function of h by 
Eq. (2-3). 


2.2.2 — Correlation Time . The RAKE receiver correlation process re- 
quires a correlation bandwidth that Is sufficiently wide to pass all 

expected doppler frequency shifts of the rain return signal. If equiva- 
lent low-pass. In-phase and quadrature-phase correlator outputs are em- 
ployed to reduce the data sampling rate to a minimum, then 


W - 


c 



(2-4) 


where fjg, Is the maximum doppler frequency, positive or negative, ob- 
served at the correlator outputs. A correlation time T^ Is defined as 


T 

c 



(2-5) 


The repetition period of the pseudorandom sequence Is To; therefore, the 
product WT - f^To defines the length of the sequence In terms of the 
number of clock periods or so-called "chips" Involved In generating one 
complete sequence. Similarly, WT,. Is the number of chips Involved In 
the RAKE correlation process. 


2.2.3 — Clutter Relectlon . The rain return signal from a selected al- 
titude region above the earth's surface Is recovered through the correla- 
tion process. However, the returns from other rain regions and the very 
strong return from the earth's surface, although at differing time delays, 
are not totally rejected In the correlation process and therefore appear 
at reduced levels as a part of the desired output. The degree to which 
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the tmdeslred, or clutter, slgnel powers can be reduced is dependent 
upon a clutter rejection factor that Is related to both WTo and 

WTq. Llndholm [4] develops an equation describing the second moment of 
the clutter level of the pseudorandom sequence correlation process for 
the case where WTc ^ WTo* These results were employed to develop an 
equation for Wy of the form 


WTo - WT +1 

W„ - ^ ^ — (2-6) 

T (WToHWT^) 

again for the region where WT^ ^WTo* This latter equation Indicates a 
best clutter rejection factor value of (1/WTo)^ is obtained when WT^ 

- WTo* As WT^ becomes Increasingly greater than WTo, If appears that 
W^ tends to follow an envelope somewhat similar In form to sln^x/x‘ 
with minimum values of (l/WTo)‘ rather than zero, asymptotically ap- 
proaching a constant value of (1/WTo)^* Best clutter rejection Is ob- 
tained with To ~ therefore To, the repetition period of the pseudo- 
random sequence. Is defined by Eq* (2-5)* 


The pseudorandom sequence parameters f^ and To as defined by Eqs. 

(2-3) snd (2-5), respectively, can be combined to determine the sequence 
length as 

■ iElfj 

and the best clutter rejection factor W_ as 

iin 

w - ^ 

^ (f To)" f" 

c c 

In order for the bilevel pseudorandom sequence to be developed as a maxi- 
mal length shift-register sequence [4], It Is necessary that 

f^To ” 2"-l, where n Is a positive Integer* It may be necessary to 

slightly adjust either f(* or To, preferably f,., to satisfy this lat- 
ter equality. It should also be noted that To must be maintained at a 
greater value than the maximum delay range of the signal returns (on the 
order of 200 microseconds) to prevent time ambiguity problems. 


(2-7) 


(2-8) 


2.2.4 — Transmission Duty Cycle . For a satellite RAKE system Involv- 
ing a single satellite and a common transmlt-recelve antenna. It will be 
necessary to alternately gate the transmit and receive functions. The 
maximum usable transmit period Is equal to the satelllte-to-earth round 
trip delay T^ where 


T 

r 


2r 

c 


(2-9) 
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which Is approximately 250 milliseconds. A SOZ duty cycle will provide 
a receive period or measured Interval of length Tf Spectral processing 
resolution will be limited to 1/T^ corresponding. In this case, to about 
4 Hz. For the satellite RAKE system, can always be greater than the 

correlation time T^.; therefore, the transmitter and receiver pseudorandom 
sequence generators can operate continuously with no dependence upon the 
transmission duty cycle. 

2.2.5 — Doppler Frequency Shift . In the previous equations of this 
section, the only parameter relating to rainfall Is fjot maximum 
doppler frequency shift. Raindrops moving toward or away from the satel- 
lite with velocity V will cause a carrier frequency shift according to 
the equation 

fd - ^ (2-10) 

where X Is the wavelength of the transmitted carrier frequency. Godard 
[5] Indicates that the maximum fall velocity Vg, of raindrops Is on the 
order of nine meters per second. Using this value for V,,, Figure A-2 
shows the relationship of 
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Figure A-2. Maximum doppler frequency shift as a function of wavelength. 
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Eccentricities of the geostationary satellite orbit can also cause a 
doppler frequency shift to all returning signals. ATS-5, for example, 
has a maximum earth radial velocity on the order of 41 meters per second 
corresponding to an 820 Hz frequency shift for a 10 cm wavelength signal 
Considering that any geostationary satellite orbit will have some eccen- 
tricity, satellite-caused doppler shifts must be expected and the satel- 
lite RAKE system Implemented accordingly. Since the satellite doppler 
shifts the earth return and rain return signals by essentially the same 
amount. It Is reasonable to assume that the effects of satellite doppler 
can be removed from the desired rain return by tracking the earth return 
signal. 


2.2.6 — Antenna Beam Width and Beam Filling Factor . The 3 dB beam 
width of a parabolic reflector antenna Is a function of both the reflec- 
tor diameter Da and signal wavelength X as shown by the following 
equation [6]: 

0 = 4, = 1.24 X 10"^ ^ (2-11) 

a 

where 6 is the antenna 3 dB beam width expressed In radians, with X 
and expressed In centimeters and meters, respectively. 

The diameter S In kilometers of the area of the earth's surface il- 
luminated by the antenna beam width Is 

S = 2r tan 6/2 (2-12) 


where r Is the satelllte-to-earth range In kilometers. With r equal 
to 3.59 X lO"* km, Eqs. (2-11) and (2-12) result In 

S = 4.45 X 10^ ^ . (2-13) 

a 

The beam-filling factor Is dependent upon the relative areas of the 
rainfall region and the antenna beam projection at or near the earth's 
surface. The area of the rainfall region of concern (7] will typically 
have diameters In the range of 20 to 50 km. If the rainfall region has 
an effective diameter In kilometers of then, with Eq. (2-13), the 

beam-filling factor Is defined as 




5.05 X 10 


r a 


D 

r 


< S 


ip = 1; > S . (2-14) 

In defining ip through Eq. (2-14), It has been assumed that some Inde- 
pendent means Is available to align the antenna such that the rainfall 
region is fully included within the antenna beam. Figures A-3a, 3b, 3c, 
respectively. Illustrate the antenna angular beam width 6, antenna 
earth coverage diameter S, and beam-filling factor ij), as functions of 
wavelength and antenna reflector diameter. It should be emphasized that 
tp cannot exceed unity. 
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Figure A-3. Antenna beam width and beam-filling factor relationships. 




2.2.7 — Atmospheric Attenuation . The equations for 

(Eqs. (2-1) and (2-2)), each contain an expression for 
atmospheric attenuation of the form 



and 


P 


rE 


A 


a 


10-0.2 jr'kdr 


(2-15) 


where k, expressed In dimensions of dB/km, is some function of the path 
characteristics and wavelength with r being the path length. Three 
sources of attenuation are significant In the wavelength range being con- 
sidered (1 cm < A ^ 10 cm). These sources are as follows: 

a. Attenuation by atmospheric oxygen and water vapor absorption. 

b. Attenuation In clouds by water droplets less than 100 microns In 
diameter. 

c. Attenuation due to scattering by raindrops. 

The integral expression of Eq. (2-15) can be represented as a sum of Inte- 
grals for each source of attenuation. In equation form. 


r r r r 

/ kdr = fkdr+ / kdr+f kdr 
■’o 'o g ■’o c ■’o P 


= a + a + a 

8 c p 


(2-16) 


where k , k , and k_ are the attenuation rate factors for gases, 
g» c’ p 

clouds and precipitation, respectively, and a , oi , and a represent 
their respective attenuation values In dB. ® ^ ^ 


The absorption by gases Is essentially Independent of meteorological 
conditions depending only upon wavelength and path length. The curve of 
Figure A-4, labeled shows the attenuation of atmospheric gases as a 

function of wavelength for a vertical path through the atmosphere. These 
values were computed for a uniform atmosphere 8 km thick with uniform 
characteristics throughout. These values closely correspond to values 
given by Hogg [8] although calculated on different bases. 

Attenuation within clouds Is caused by water droplets with radii of 
100 microns or less. The amount of attenuation Is dependent upon the 
liquid water content of the cloud. The curve of Figure A-A, labeled a^,, 
was obtained using data from Bean and Dutton [9], p. 291. The assumed 
meteorological conditions for this calculation are graphically illustrated 
by Figure A-5. This figure also shows the typical conditions assumed 
throughout this report except where otherwise stated. 

The attenuation due to rain Is a function of rainfall rate. The 
values for rain attenuation shown in Figure A-4 as Op are from data for 
an average rainfall rate of 25 mm/h [9], and assumes a 1 km thickness of 
rain above the observation region. 
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o,c CLOUD 



Figure A-4. One-way signal attenuation 






*R USED FOR CALCULATION OF 
RAIN ATTENUATION ONLY 


Figure A-5. Storm cloud model for signal attenuation calculations. 
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RELATIVE POWER LEVEL 
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Figure A-6. Relative power return versus raindrop terminal velocity. 
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2.2.8 — Rain Doppler Spectrum Estimate . In order to determine a 
value for minimum detectable received signal power, it is 

necessary to have some estimate of the doppler spectrum of the rain re- 
turn signal. Each raindrop will have a radar cross-section a propor- 
tional to D®, where la the drop diameter. The total reflected 

power on a relative basis can then be expressed as 

N 

d‘ (2-17) 

^ 1=1 

where N Is the total number of drops In the observation region. By 
grouping the drops Into Intervals according to their diameter, Eq. (2-17) 
can be written as 


m P N 

’’r “ ( 2 - 18 ) 

where Pj Is the percentage of the total number of drops that fall Into 
the jth Interval of diameter and m Is the number of intervals. 

Figure A-6 Illustrates the relative rain-power distribution as a func- 
tion of raindrops terminal fall velocity. This figure Is a plot of Eq. 
(2-18) with Pj data obtained from Table 1. Ci was chosen to give powers 
groportlonal to R^'® (R = rainfall rate In mm/h) as in the equation for 

P . In addition, drop velocities were derived from drop diameter using 
raR 

the relationship [10] 


Vd = 13v^ 


(2-19) 


It Is apparent from Figure A-6 that the half-power response points corre- 
spond to a three- to four-meters per second spread, regardless of the 
rainfall rate. Designating this velocity spread value as AV, the 
equivalent doppler frequency spread (or doppler bandwidth) can be 
obtained from 


B 


d 


2AV 

A 


( 2 - 20 ) 


Assuming the four meters per second velocity spread value. 


B 


d 


800 

A 


with A expressed In centimeters. 


( 2 - 21 ) 


2.2.9 — RAKE Receiver Sensitivity . The minimum detectable received 
power Is a function of both the effective noise temperature of the sys- 
tem receiver, Tsys, and the rain doppler spread or bandwidth Bj, where 
Bj Is defined by Eqs. (2-20) and (2-21), and 
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Table 1 


Raindrop Size Distributions* 

(Percent of volume contained drops of size D^) 


Drop Size 



Rainfall Rate 

R (mm/hr) 


(D, mm) 
a 


2.5 

25 


100 

.5 


7.3 

1.7 


1.0 

1.0 


27.8 

17.6 


4.6 

1.5 


32.8 

18.4 


7.6 

2.0 


19.0 

23.9 


11.7 

2.5 


7.9 

19.9 


13.9 

3.0 


3.3 

12.8 


17.7 

3.5 


1.1 

8.2 


16.4 

4.0 


0.6 

3.5 


11.9 

4.5 


0.2 

2.1 


7.7 

5.0 



1.1 


3.6 

5.5 



0.5 


2.2 

6.0 



0.2 


1.2 

6.5 





1.0 

7.0 





0.3 

* 

From Radio Meteorology, 

Bean and Dutton NBS Monograph 92, 1966, 

p. 296. 
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sys A L 

(1-a) + T + ;[ 
m=2 

Vi 

(2-22) 

with 






a 

= 

Transmission line loss between antenna and 

receiver 


= 

Effective antenna temperature 




= 

Transmission line temperature 




= 

Effective noise 

temperature of the first stage amplifier 

T 

m 

a 

Effective noise 

temperature of the mth stage amplifier 

Vi 

= 

Gain of the (m-l)st stage amplifier. 



For the satellite RAKE application, and are both reasonably 

assumed equal to 290°K and only the first two amplifier stages are con- 
sidered such that 

T = 290 + T, + . (2-23) 

sys ^ ^1 

The minimum detectable received power P^j, Is assumed to be equal to the 
noise power of the receiver system within the detection bandwidth. The 
receiver noise power Is determined by the equation 
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(2-24) 


P 

n 


k T if 


sys 


where 


—23 

k = Boltzmann's constant = 1.38 x 10 joules/°K 
"*■878 "" Effective receiver noise temperature (°K) 

Af >= Noise bandwidth (Hz) 


Curves of Prm “ Ejj are shown In Figure A-7 for various bandwidths, 

B = Af as a function of wavelength or frequency, considering the noise 
temperatures of receiver front-end devices tht are practical for satel- 
lite use.* The curve of Figure A-7 for B = Bj represents the minimum 


detectable rain return power, • Eor the frequency range from 3 to 

18 GHz, field-effect transistor (FET) or tunnel diode preamplifiers are 
employed. Above 18 GHz, mixer type front-ends are used. More sophisti- 
cated preamplifiers, such as parametric amplifiers and masers, can pro- 
vide lower receiver noise temperatures; however, these devices are not 
considered suitable for satellite applications at this time. 


2.2.10 — Empirical Correction Factor . The correction factor F used 
In the basic rain return equation Is based on data from Battan [1], p. 60. 

A linear relationship with X was assumed In the plot shown in Figure A-8. 
This correction factor Is purely an empirical relationship and is based on 
the discrepancies between observed and calculated rain return powers. 


2.2.11 — Reflectivity Factor and Rainfall Rate . The reflectivity fac- 
tor Z Is equal to the sum of the reflectivity of each drop throughout 
the rainfall volume, as discussed In paragraph 2.2.8. Data have been col- 
lected on the relationship between Z and the rainfall rate R. Battan 
[1], p. 56, gives 23 different relationships for various types of rainfall 
and locations. The representative equation that best fits these data and 
that will be used throughout this study Is: 

Z = 2 X 10"^° R^-^° (2-25) 

where Z has dimensions of cm*/cm^ with R in mm/hr. 


2.3 Optimization Procedures 

Two criteria for optimization were considered. It was first assumed 
that the system performance would be limited by the presence of an earth 
return signal of sufficient strength to obscure the rain return signal. 
Earth return-to-raln return power ratios of up to the order of 90 dB can 
be expected. Assuming the RAKE correlation process can reduce the power 
level of the earth return by at least 40 dB, the first criterion was based 


Receiver noise figures ranging In value from 4.5 dB at X = 10 cm to 
12 dB at X = 1 cm were assumed. These variations with X are reflected 
in the graphs In Figure A-7. 
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Figure A-8, Correction Factor F versus X* 


upon the capabilities of the receiver (after correlation) and spectral 
processing to separate the desired rain return from the reduced earth re- 
turn for a rain return-to-reduced earth return power ratio of -50 dB 

or 10 “ 5, 

The earth return power will be reduced by the clutter rejection 

factor Wx; therefore the desired rain return-to-reduced earth power ratio 
will be 


raR 


(2-26) 


Starting with basic Eqs. (2-1) and (2-2) for and respectively, 

and using Eq. (2-8) for and appropriate substitutions of Eqs. 

(2-7), (2-10), (2-14) and (2-25), an equation defining the minimum de- 
tectable rainfall rate for the constraints imposed by Eq. (2-26) can be 
derived as 
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min 


R, [cot^Soexp (-tan^a/tan^So) ] 


10.625 


hX' 


(2-27) 


(8 


.76 X 10"^)lt®F|K|^10“°'^/j ’’r 


2 

a 


This latter equation is plotted in Figure A-9(a) as a function of X for 

= 10~^, 8o =5°, a = 0°, |k|^ = 0.92, h = XO' cm, = 50 km and 

Dg « 10 in with the atmospheric attenuation and F factors as shown in 
Figures A-4 and A-8, respectively. 

Now, by substituting R^in ^*1* (2-27) into the rain return Eq. 
(2-1) and through use of the relationships defined by Eqs. (2-11) and 
(2-25) plus the previously undefined relationship between antenna aper- 
tural area and antenna reflector diameter. 


A = 7.85 X 10^ 

P a 


(2-28) 


an equation defining the minimum transmitter power required to observe 
the previously related minimum detectable rainfall rates is derived as 


.10 


.0.2L,s 


raR 


min 


*'mln D* R^h^[cot^6o exp(-tan*oc/tan^8o) ] 


(2-29) 


''hefe (PraR> 


Is the minimum detectable rain return power in band- 


mln 


width Bj shown in Figure A-7. Equation (2-29) is shown plotted in 

Figure A-9(b) as a function of X for antenna diameters of 10 m and 20 m, 
L = 2.5 dB and r = 3.59 x 10® cm with the remaining parameter values as 
employed in determining 

The second criterion for optimization was simply to determine the 

transmitter power required so that the rain return power, P , would 

raR 

exceed the receiver threshold sensitivity, (P ) , for a constant 


raR 


min 


value of rainfall rate R. In this case, starting with Eq. (2-1) for 

p , substituting the relationships defined by Eqs. (2-11), (2-14), 
raR 


(2-25), and (2-28), letting P 


raR 


> 

min 


and solving for P^, 


0.483 A'r-io”-^^^,^^^)^^^ 


TT® D‘*D'‘hF|KPR^'*10 




(2-30) 
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Pt from Eq. (2-30) Is plotted In Figure A-10 as a function of \ for 
rainfall rates of 10, 30 and 100 mm/hr with r ° 3.59 x 10^ cm, ° 10m, 

Dj = 50 km, (k|* = 0.92, h = 10® cm and L «■ 2.5 dB. Attenuation fac- 
tor (P ) and F are, respectively, as defined In Figures A-4, A-7 
rsR min 

and A-8. 


2.4 Discussion of Optimization Results 

The curve of Figure A-9(a) Indicates that for a rain return-to-reduced 
earth return ratio of 10~^ and with the other assumed conditions, an ab- 
solute detectable rainfall rate theoretically occurs at a wavelength of 
approximately 2 cm, or 15 GHz. However, to observe the slightly 

more than 0.2 mm/hr value of rainfall rate at that wavelength requires a 
transmitted power of more than three kilowatts when using a 20-meter 
diameter antenna reflector. Higher values of rainfall rate can be ob- 
served using lower transmitted powers. However, lower rainfall rates can- 
not be observed If higher transmitted powers are employed because of the 
rain return-to-reduced earth return ratio limitation. One obvious way to 
reduce the transmitted power requirement would be to Increase the antenna 
diameter beyond 20 meters. Antenna size cannot be Increased without limit, 
however, because the earth coverage area decreases as antenna size In- 
creases. Too small a coverage area Increases the problems of obtaining 
Integrated rainfall rate data. Also, since antenna diameter really Im- 
plies antenna gain there Is a practical limit to antenna gain due to the 
precision with which the parabolic reflector can be constructed. A 
further discussion of practical limitations Is presented In paragraph 2.5. 

The curves of Figure A-10 Indicate that optimization on the basis of 
receiver sensitivity Is the more realistic approach. Again, the minimum 
transmitted power required Is theoretically obtained at a wavelength of 
about 2 cm regardless of rainfall rate If the rainfall rate Is above the 
minimum Indicated In Figure A-9. In order to observe a rainfall rate of 
10 mm/hr using a 10-meter antenna, a transmitted power at 2 cm of 100 
watts is required. For a 20-meter antenna, the transmitted power Is 
reduced to about 7 watts. 


2.5 Practical Limitations 

In the previous discussions the values of antenna gain, as related to 
antenna diameter, and transmitter power have been allowed to vary essen- 
tially without limit. The results obtained on this basis tend to Indi- 
cate the geostationary satellite RAKE concept for measuring rainfall rate 
Is theoretically feasible, assuming a reasonable Interpretation of the 
rainfall doppler spectrum Is possible. Before the concept can be con- 
sidered feasible on a more practical basis, however. It Is necessary to 
consider the effects of antenna gain and transmitted power limitations 
Imposed by present and near future state-of-the-art technology. 

The satellite antenna gain that can be obtained Is limited by the 
precision with which the parabolic reflector can be constructed and as- 
sembled In space. Earth antennas have been constructed that have 70 dB 
and more of gain In the 3 to 15 GHz range. However, the antennas for the 
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Figure A-10. Transmitter power required versus X, R. 
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lower frequencies are several hundred feet In diameter and at the higher 
frequencies are precision ground solid reflectors. It Is known that 
satellite antennas with gains up to 54 dB are under consideration for some 
future space stations, but these are on the order of IS ft In diameter 
and are not folded for launch. For the present and Immediate future a 
50 dB antenna gain for the 3 to 15 GHz range seems to be a reasonable 
limit. In the near future, antenna gains of 60 dB may be possible In 
satellite applications. It Is also considered chat for both the pres- 
ent and near future, satellite antenna diameters will be limited to about 
20 meters, particularly for wavelengths In Che 2 to 5 cm range. 

The most desirable method of power amplification Is solid state be- 
cause of Its ruggedness, high reliability, low weight and high efficiency. 
However, the present state of microwave semiconductors limits the power 
generation to about 200 watts at 1 GHz and 402 efficiency, declining to 
about 30 watts at 4 GHz and 102 efficiency. 

The second choice for power generation would be traveling wave tubes 
(TUT's) because of their moderate ruggedness and weight. However, their 
efficiency Is only 20-252 and the gain Is low requiring more preampllf Ica- 
tlon. TVJT's that generate 10 kW at 1 GHz and up to 1. kW at 10 GHz are 
available. However, any over several hundred watts require water cooling 
which becomes Impractical for anything but a very large satellite. 

Very high powers, 100 kW and above, can be generated by klystrons In 
the 1 to 10 GHz range. However, the klystron Is the least rugged of the 
devices and has about 302 efficiency. As with the TKT, water cooling Is 
required for powers above several hundred watts. 

For the present, a TWT limited to about 250 watts of power appears to 
be the best choice for a transmitter power amplifier. As satellite avail- 
able primary power and weight capacities Increase, transmitter output 
power can also Increase. In the near future, transmitter powers on the 
order of 2.5 kW should be possible. 


To best Illustrate the effects of the limitations on antenna gain and 
transmitter power as discussed above, the following equation defining 
rainfall rate was developed: 


(2.56 X 10^)A^r‘*10 


0.2(L + / kdr) 

° (p ) 

' raK'mln 


Tr’p^G^D^hF|K|^ 


'|0.625 

► 


(2-31) 


where G Is the numerical antenna power gain and all other parameters 
are as previously defined. Equation (2-31) has been employed to deter- 
mine the two curves of rainfall rate versus wavelength Illustrated by 
Figure A-11. These two curves show what Is considered to be the pres- 
ent and near future capabilities of a practical satellite RAKE system. 


Figure A-11 Indicates that the current state-of-the-art does not per- 
mit detection of rainfall rates below about 150 mm/hr. This Is not 
considered suitable. However, with anticipated Improvements In antenna 
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gain and transmitted power capabilities, measurement of rainfall rates In 
the more desirable 10 to 25 mm/hr range will become possible. Both curves 
of Figure A-11 Indicate wavelengths between 2 and 6 cm offer the greatest 
capability. 


3. Techniques for Implementation - 

The following discussion covers methods of hardware implementation of 
the rake system to obtain the theoretical performance Indicated In Sec- 
tion 2. It will consider system problems In the Implementation of the 
hardware necessary to obtain the data, and will not dwell on circuit de- 
sign. Three approaches to the data collection will be considered. In 
the first, all of the RAKE equipment is in the satellite with the correla- 
tion data telemetered to a ground control station. In the second, the 
RAKE equipment Is In a ground station with only a transponder in the sat- 
ellite. The third approach uses separate satellites for the RAKE trans- 
mitter and receiver to obtain antenna Isolation and thereby eliminate the 
need for transmlt-recelve signal gating. 


3.1 All Measurement Equipment in the Satellite 

The first approach to be considered is to place all of the data col- 
lection equipment In the satellite. Commands from the ground would be 
used to Initiate the experiment and aim the antenna at the desired ground 
location. The detected return from the rain would be telemetered back to 
ground station for further analysis. 

Figure A-12 Is a block diagram of the RAKE hardware that would be In 
the satellite. A single parabolic antenna would be used by both the RAKE 
transmitter and receiver. The use of separate antennas to transmit and 
receive does not. seem practical when considering a satellite perhaps the 
size of ATS-F. It Is of more Importance to share as large an antenna as 
practical rather than have separate but smaller antennas. Since it Is 
not possible to obtain sufficient Isolation between a transmitter and 
receiver sharing the same frequency band, the transmitter must be turned 
off while receiving. The two-way propagation delay from the satellite to 
earth Is approximately 250 milliseconds. Thus a cycle of transmitting 
for 250 milliseconds and receiving for 250 milliseconds would be reason- 
able. Shorter cycles could be used but they do not offer any advantages. 
Depending upon frequency, some form of circulator would be used to 
direct the transmitter power to the antenna and the receive power to the 
receiver. 

The transmitter frequency would be a multiple of a stable crystal os- 
cillator with the pseudorandom sequence modulated on the signal. The 
degree of modulation would be such that after multiplication the modula- 
tion would be + 90 degrees. If the pseudorandom sequence were perfectly 
balanced, there would be no carrier remaining for exactly + 90 degrees 
modulation. 

Since the satellite can have a vertical component of velocity there 
will be a doppler on both the return from the earth and from the rain due 
to the satellite motion. To remove this doppler from the measurement, 
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Figure A-12. Satellite RAKE transnitter/recelver , block diagram. 

















the receiver local oscillator la phase- locked to the earth return. As 
shown for the earth return channel In Figure A-12, the received signal Is 
mixed with the local oscillator which has been modulated with a replica of 
the transmitted pseudorandom sequence. When the two sequences are In 
synchronism, the carrier Is regenerated and amplified. The phase detec- 
tor and VCO complete the loop and maintain the local oscillator In syn- 
chronism with the earth return regardless of satellite velocity. 

The vertical velocity of the satellite also causes a doppler shift In 
the clock frequency of the received pseudorandom sequence. Compensation 
must be made for this received clock frequency shift in order for the se- 
quence generated at the receiver to stay In synchronism with the received 
sequence. This can be done by making the clock frequency a fraction of 
the carrier frequency. If at the transmitter the clock is 1/Nth of the 
carrier, the doppler on the clock in the received signal will be 1/Nth of 
the doppler on the carrier of the received signal. When the locally gen- 
erated carrier at the receiver Is phase-locked to the received carrier, 
the carrier can be divided by N to obtain a clock for the receiver se- 
quence which Is In phase with the doppler-shlf ted received clock. This 
technique for generating sequence clocks is shown in Figure A-12. 

The design of the rain return channel presents special problems caused 
by the large earth return signal. Paragraph 2.3 of this report has Indi- 
cated that the earth return power can be 90 dB greater than that of the 
rain return. Any portion of the receiver handling both signals simul- 
taneously must be linear over the 90 dB range to prevent Intermodulation 
and harmonic distortion. Otherwise the ability to select a small rain 
signal out of the clutter will be greatly reduced Impairing the success 
of the measurement. 


For this reason the correlation of the return signal with the local 
pseudorandom sequence Is performed In the first mixer before any amplica- 
tion to obtain some rejection of the earth return immediately. The mixer 
can be designated to be linear over this range of signals as long as the 
total signal power Is not greater than about -30 dBm. 


The bandwidth of the IF amplifier after the mixer must only be wide 
enough to pass the doppler spectrum of the rain return which Is within 
+ fjjjjj from the center frequency of the return. The doppler fj Is 
caused only by the rain since the satellite doppler has been removed by 
phase-locking on the earth return. The spectrum of the signal in the IF 
Is anticipated to be essentially as shown In Figure A-13. The earth re- 
turn has been reduced to a residual carrier and line spectrum, while the 
rain return Is offset by the rain doppler. Section 2 has shown that the 
residual carrier portion of the earth return can be rejected by 


)‘, or about 40 dB, by the correlation process. After selecting only 

the return at + f<j, the earth return is still about 50 dB above the rain 
return. A rejection filter in the center of the IF band can further re- 
duce the earth return so that the linearity requirement in the remaining 
portion of the receiver Is minimal. 


Two phase detectors In quadrature are used for the final detection of 
the rain return. The final filters are each fj wide to pass the rain 
but ac coupled to reject the remaining earth return at dc. 
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Figure A-13. IF spectrum. 


The same pseudorandom sequence generator Is used In the earth and 
rain channels but Is delayed for the earth channel by the propagation time 
corresponding to the altitude above the earth of the rain segment being 
examined. Once the earth loop Is phase-locked, the altitude of the rain 
segment examined can be selected by changing the relative delay of the 
sequence applied to the rain channel. This can be controlled by the com- 
mand link from the ground. 

The problem of acquisition of the earth return requires two steps: 
first, locking on the carrier and then onthe pseudorandom sequence. Suf- 
ficient carrier power must Initially be left In the transmitted signal to 
obtain phase lock. Then the sequence can be stepped until sync Is ob- 
tained. After complete acquisition, the residual carrier power In the 
transmitter can be reduced. 

Table 2 shows the expected losses for a telemetry link from the sat- 
ellite to the ground station. It was assumed that a 30-foot steerable 
parabolic antenna on the ground can be used for the link. For Che lower 
frequencies considered In Che table, little antenna gain can be obtained 
in the satellite antenna without using a very large array or a parabolic 
reflector. Thus a dipole antenna was considered appropriate for the 
lower frequencies. Since It Is linearly polarized it will be necessary 
to use a circularly polarized ground antenna. This results In a 3 dB 
additional loss. 

For the midfrequencies a helical antenna can be used. This antenna 
provides Increased gain over Che dipole and eliminates the polarization 
loss. A helical antenna length of 3X/2 is typical of the size. 

For the higher frequencies, a horn antenna can be used. For this 
application the gain of the horn antenna is limited so that the 3 dB 
beamwldth will be on the horizon even when the satellite Is tilted to 
point to the opposite horizon. The horn antenna Is also linearly polar- 
ized. 
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Table 2 


Transmission Loss — Telemetry Signal 


Freq 

MHz 

Isotropic 

loss 

(dB) 

30-ft 

ground 

antenna 

gain 

(dB) 

Satellite 

antenna 

gain 

(dB) 

Polarlza 

tlon 

loss 

(dB) 

Beam 

edge 

loss 

(dB) 

Total 

trans- 

mission 

loss 

(dB) 

100 

165 

17 

Dipole - 

2 3 

_ 

149 

200 

171 

23 

Dipole - 

2 3 

- 

149 

500 

179 

31 

Helical - 

10 - 

1 

139 

1,000 

185 

37 

Helical - 

10 - 

1 

139 

2,000 

191 

43 

Helical - 

10 - 

1 

139 

5,000 

199 

51 

Parabola 

- 14 3 

3 

140 

10,000 

205 

57 

Parabola 

- 14 3 

3 

140 


Table 3 

Telemetry Power Requirement 

^d 

(Hz) 

40 X f 
(dB) 

Trans 

loss 

(dB) 

S/N 

(dB) 

P„ PER Hz 
N 

(dBW) 

"t 

(dBW) 

100 

36 

-139 

+10 

-200 

-15 

300 

41 

-139 

+10 

-200 

-10 

1000 

46 

-139 

+10 

-200 

- 5 


Table 2 shows that the selection of telemetry frequency Is not criti- 
cal above 500 MHz. Thus the selection can be made based upon other fac- 
tors such as the availability of equipment. 

The bit rate for the data will depend upon the doppLer frequency. 
Assuming a sample rate of twice the doppler (Nyqulst rate) and an analog- 
to-dlgltal converter of 10 bits, the two channels require 40 fj bits per 
second . 

The data can best be transmitted using coherent phase modulation. 

For a 10”^ bit error rate, a 10 dB slgnal-to-nolse ratio In a bandwidth 
of 1 bit rate Is required. Assuming a 7 dB receiver noise figure. Table 
3 shows the required transmitter power. The values are small enough so 
that no problems should be encountered In supplying the required power. 
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3.2 Satellite Transponder 

The second approach to be considered Is to have the RAKE transmitter 
and receiver on the ground and a transponder In the satellite. The sat- 
ellite equipment for this approach Is shown plctorlally In Figure A- 14. 

While this approach greatly simplifies the satellite hardware, It has 
one problem that disqualifies It from further consideration. This Is the 
fact that the earth return Is on the order of 90 dB above the desired 
rain return, and linear amplification Is required until correlation Is 
performed. Obtaining a hlgh-power amplifier with this characteristic Is 
not feasible. In addition, the satellite power for relaying the signal 
to the ground station Is beyond practical limits. For example, assuming 
a 300 Hz rain doppler and a 7 dB noise figure, the required satellite 
power Is: 


Holse power per Hz 

= -200 

dBW 

Bandwidth (S “ N) 

= 25 

dB 

Required signal power 

-175 

dBW 

Loss (Table 2) 

= -139 

dB 

Trans signal power 

= - 36 

dBW 

Earth return to signal ratio 

90 

dB 

Total transceiver power 

+ 54 

dBW, 


or 600,000 watts 


This approach will not be considered further. 


3.3 Two-Satellite Configuration 

The third approach to be considered Is to have the RAKE transmitter 
In one satellite and the receiver In a second satellite. The division of 
equipment Is obvious from Figure A-12. The advantage of this arrangement 
over the single satellite Is that transmission can be continuous rather 
than have a 50% duty cycle. Continuous transmission will allow greater 
resolution In the frequency analysis of the rain return. For example, a 
sample 250 milliseconds long limits resolution to approximately 4 Hz. 
However, the cost of obtaining the Improved resolution Is almost double 
the single satellite approach and does not appear to be justifiable. 


4. Implementation In Planned Geostationary Satellites 

In reviewing planned satellites In which the RAKE system might be 
Implemented, the primary characteristics required were as follows: 

(a) Geostationary orbit, 

(b) Low orbit eccentricity, 

(c) Antenna with a gain of 50 dB or better at 3 cm wavelength. 
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(d) Attitude stabilized so antenna pointing of about + 0,15 degree 
(1 dB width) can be maintained, 

(e) Space available for the RAKE system. 

Inquiries were made to Goddard Space Flight Center, Comsat Corpora- 
tion, and several private contractors to see what satellites presently 
being planned would be available and would come close to meeting the pri- 
mary requirements. The only nonmilitary satellites approaching the re- 
quirements appear to be ATS-F and ATS-G. To our knowledge the design of 
ATS-F which Is to be launched In 1974 has been frozen and Is no longer 
open to new experiments. ATS-G design Is still open, but some persons 
have expressed doubt that It would ever be built. 

The antenna planned for ATS-G Is a 10 meter parabola with a peak gain 
of about 49 dB at 6 GHz, and a decreasing gain above this frequency. The 
peak available primary power Is about 450 watts so the transmitted power 
obtainable would be about 125 watts. Thus, If ATS-G were used, the values 
of present state-of-the-art rainfall detection shown In Figure A-11 would 
have to be Increased slightly to conform with the power and antenna gain 
available. With a transmitter frequency of 6 GHz and a transmitter power 
of 125 watts, the minimum detectable rainfall Is about 225 mm/hr. 

The ATS-G satellite has a ground-to-satelllte control system for atti- 
tude control and experiment control that appears to be adequate for the 
rake program. It also has a telemetry communication system that would be 
adequate. The RAKE transmitter, receiver and correlation equipment would 
have to be developed especially for the satellite as their present trans- 
mitters and receivers are not adequate. 


5 . Conclusions 


The results of the system analysis and optimization Indicate that the 
measurement of rain doppler spectrum Is feasible using the geostationary 
satellite RAKE concept. A major problem does appear to exist, not in ob- 
taining the rain doppler spectrum, but rather In Interpretation of this 
spectrum In terms of rainfall rate. Figure A-6 offers a clear Illustra- 
tion of the potential problem. In this figure the doppler spectrum sig- 
natures for the three rainfall rates differ little In either frequency 
value or frequency spread. They do differ more widely in level. How- 
ever, level variations are related to both rainfall rate and beam-filling 
factor. It was considered necessary to Indicate the presence of this po- 
tential problem, although no solution to the problem was pursued during 
the course of this study. 

State-of-the-art limitations to satellite antenna gains and trans- 
mitter powers do limit the practicality of employing the geostationary 
satellite RAKE concept to obtain rain doppler spectra at the present time. 
The future shows good promise for the concept, however, as satellite 
antenna gains and transmitter powers are increased. 

The portion of the study pertaining to methods of implementing the 
satellite RAKE system has resulted in a clearly defined "best" method. 

This method of Implementation requires that all RAKE equipment be Incor- 
porated Into a single satellite with only control signals and RAKE 
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correlator output data linked to the ground based terminal. ATS-G appears 
to be the only nonmilitary satellite to be placed In orbit at some near 
future date that has capabilities approaching those required for satellite 
RAKE Implementation. However, an ATS-G Implementation would not be suit- 
able unless greater transmitter powers can be made available. 

If the problem of Interpreting rain doppler spectra in terms of rain- 
fall rate can be circumvented, the future capabilities of the satellite 
RAKE system to measure rainfall rate would appear to be assured. 
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POINTING ERROR ANALYSIS OF GEOSYNCHRONOUS SATELLITES 


Anlruddha Das and T. C. Huang 


FOREWORD 


Introduction 


The two principal problems In unmanned satellite dynamics are: 

1) stability of the satellite; and 2) precise determination of the orienta- 
tion of the different sensors mounted on the satellite. In the following 
analysis the latter problem Is treated. 

To recognize the problem, we have to note the following facts. A sat- 
ellite usually has two sets of sensors. The first set, the Star Tracking 
Sensors (STS) , looks at some preassigned stars and sends down the signal 
for the angular positions of the stars relative to the sensors. From these 
data, the attitude angles of the satellite are calculated. The second set 
looks at the Earth and Is called the Earth Viewing Modules (EVM) . For 
scientific purposes. It Is obviously necessary to know the location of the 
point on the Earth at which the EVM was looking when the data were sensed. 

So far, the practice has been to calculate the attitude angles of the EVM 
from the attitude angles of the STS. This would have been satisfactory 
If both sets of sensors were part of the same rigid body, but In reality 
the flexible elements on which these sensors are mounted lead to errors 
In the EVM attitude angle determination. 

Another problem Is that the data are sensed by both the STS and EVM at 
discrete time Intervals. Also, the sensing time and frequency are not 
necessarily synchronous for the two sets of sensors, so for correct corre- 
lation of data It Is necessary to know the motion of the EVM during the 
Interval between two sets of data transmitted by the STS. To do this, an 
accurate dynamic model of the satellite Is required which takes Into ac- 
count the flexibilities, energy dissipative sources and environmental torques 
Influencing the motion. 

Now, trying to solve the dynamic equations which Involve the angular 
velocities of the satellite u leads to further complications. To obtain 
these angular velocities as a function of time In an Interval, say 0 t ^ T, 
the magnitudes of u have to be known at one point In the Interval, say at 
t = 0. These Initial values of w have to be obtained from the data sent 
by the STS. These data are Inaccurate, however, as the exact STS angular 
positions are not known. To calculate the displacement of the STS, u and 
the Initial values of the coordinates of the STS qo have to be known. 

Since there Is no way to measure qo, one can only obtain the most probable 
values of bi from a set of assumed values of qo to start with and 
then use Improved values of qo In successive Iteration processes. 
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The environmental and controlling torques are associated with torque 
noises. Another source of error Is the presence of electrical noise In the 
data transmitted by the sensors. The presence of these noises makes every 
variable In the problem a random variable. 

The solution of the dynamic problem requires the solution of the asso- 
ciated control problem. In addition, an optimal control system has to be 
formulated to meet the objectives of accurate prediction of the attitude 
angles of the satellite. 


Literature Review 


Most of the early work In satellite dynamics was done to analyze the 
stability of the motion, and It has thrown a considerable amount of light 
on this problem. 

The "major axis" rule for the stability of freely spinning vehicles was 
first published by Fllklngton [1], Bracewell and Garlott [2], and Perkel 
[3]. Most of the existing spln-stablllzed satellites were designed with 
this criterion. Quantitative determination of the Influence of energy 
dissipation by dampers on the motion of bodies spinning about their major 
axes Is shown by Thomson and Reiter [4], and Likins [5], [6]. The dual- 
spin concept was first developed by A. J. lorlllo [7], who also considered 
simple damper models In both bodies. They were closely followed by Karymov 
[8], Rossi et al. [9], Likins [10], Likins and Mlngorl [11], Mlngorl [12], 
and Pringle [13]. 

The usual technique used by these authors was to linearize the rigid 
body Euler equations and then to apply the Routh-Hurwltz criterion. Mlngorl 
[12] compared the linear predictions with Floquet theory. Likins and Mln- 
gorl [11] obtained a Liapunov function for linear systems. Pringle [13] 
derived his theorem based on Liapunov's direct method for a lumped mass 
system, after linearization. An Interesting method was shown by Flatley 
[14], [15]. He used the angular velocities from a rigid body assumption 
to solve a simple spring-mass-damper model, and obtained results very close 
to Mlngorl's Floquet analysis for a similar system. 

The most significant contributions to the problem of obtaining the mode 
shapes for a flexible satellite are from Likins and his associates [16], 
[17], [18], [19], [20]. They considered distributed mass systems by lumped 
mass approximation In [16], [18] and [19], and completely linearized the 
systems with distributed mass mode functions In [17] and [20]. The forcing 
and controlling torques were not considered. In their formulation, they 
did not obtain any equations free of the spatial coordinates. 

Their approach was fully deterministic, and the problem of nondeflnlte 
Initial conditions was not considered, f^e and Robe [21] considered two 
symmetric rigid bodies connected by a flexible beam. A rigid satellite 
with two torsion pendulums was considered by Balnxmi and others [22]. Flex- 
ible coupling of two rigid spinning bodies was also treated by Cretcher and 
Mlngorl [23] and Wenglarz [24]. In a comprehensive report prepared by 
Avco Systems Division for Goddard Space Flight Center, Maryland [25], the 
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problem of sensor noise for the angular velocities of a rigid satellite 
with spring-mass systems was considered. The dependence of the error on 
the sensor oscillations was not taken Into account. Various problems of 
control were considered In [26] for simple linear systems. 

The major area of thermally induced flutter of a general flexible satel- 
lite has not been considered thus far. The case of long beams only has 
been considered In Etkln and Hughes [27]. With the exception of [25] and 
Zach [29], the system analyses mentioned so far have left out the effect of 
environmental torques. The torques acting on a rigid satellite are con- 
sidered In [25]. Zach [29] used a simple torque model for gravity-gradient 
satellites. In addition, Dobrotln [30] and Tidwell [31] made the only 
realistic torque models, published so far, without going Into the dynamic 
system. The solution for the vibrations of a fixed base, elastic, dlstrlb- 
uted-mass structural system. Involving shells, plates, beams, rigid bodies 
and point masses, has been developed by Huang and others [32], [33], [34], 
and [35]. The present analysis, which also considers translational and 
rotational motion of such complex structures, will be an extension of that 
series of work. 

We have reviewed the literature dealing only with the general area of 
satellite dynamics. References related to the methods of analysis will be 
made when needed. 


The Scope of this Study 

In the present analysis the basic configuration of a dual-spin satel- 
lite has been made. This can be modified to a spinning or a three-axes 
stabilized satellite. The dynamic model Includes elastic shells, plates, 
beams, rigid bodies and point masses. The effects of thermal stresses, 
large angular velocities and the effect of the motion to the center of mass 
due to vehicle deformation are Included. In this formulation, the spatial 
dependences are maintained linear while the time dependences are nonlinear. 

The analysis proceeds with the following plan. We start with appro- 
priate partial differential equations for the beams, plates and shells In 
global coordinates fixed to the satellite. These are then transformed Into 
local coordinates fixed to the nominal positions of each of the elements. 
Based on the linear theory the solution spaces of the distributed mass ele- 
ments are chosen; this Is done by choosing a finite series of terms of 
known spatial functions multiplied by unknown time-dependent coefficients 
qi. Then, using Galerkln's method, the spatial functions are Integrated 
out of the partial differential equations. This leads to a matrix of 
ordinary linear second order differential equations In terms of the unknown 
coefficients of the spatial functions. The coefficient matrices of all but 
the second derivatives In these equations remain nonlinear functions of the 
angular velocities of the satellite. The angular momentum equations 

for the whole satellite are then obtained, which Involve both the sets 
q£ and . From all these equations, the unknowns qjL are eliminated 
leading to nonlinear equations In only. These equations turn out to 

be singularly perturbed equations and are scaled out to obtain short-time 
and long-time equations for the angular velocities It Is found that 
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the long-time equations, on being drastically simplified, lead to the 
usual rigid body equations of motion of satellites. The short-time equa- 
tions lead to hitherto unknown equations for the transient response of 
flexible satellites. Solution of this equation requires complete informa- 
tion of the controlling torques for the satellite, so a time-optimal control 
policy Is then assumed and a suitable optimum control system Is then formu- 
lated. A linearized solution of the system Is then subjected to statisti- 
cal analysis for getting the probability distribution functions of the 
required variables. The mean values then are to be obtained from the non- 
linear equations by a numerical procedure. 

The primary advantage of the method outlined is that Irrespective of 
the number of flexible elements In the dynamic model, we will always have 
to solve only three nonlinear ordinary differential equations In three un- 
known variables. Also the boundary conditions for the q^ are satisfied 
before solving the problems by suitably reducing the number of Independent 
elements In the set {q^}. 


ObJ ectlves 


The objective of this analysis Is to estimate the pointing error of 
different satellites. For a three-axes stabilized satellite, a very pre- 
cise determination of Instantaneous attitude angles Is possible. However, 
Its motion at a subsequent time and the corresponding control Is very un- 
certain. A spinning satellite provides a very stable platform in space, 
but the attitude measurement Is uncertain because of the structural flexi- 
bility. 

In this study one or more of the rigid bodies will be the model of 
the attitude determination sensors. The rest of the rigid bodies will be 
the EVM's. Then, for a particular design, this analysis will provide: 

a) the extremes of the attitude error between the different sensors; 

b) a probabilistic time history of the magnitude of the error for 
the transient zone after every control torque pulse; 

c) a computer program to plot out the pitch, roll and yaw limit 
cycles for a 90Z probability density bandwidth; 

d) an estimate of the stiffness requirement of the flexible elements 
for a given maximum error limit; and 

e) comparison of all these for spinning, nonspinning and dual-spin 
satellites. 
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nomenclature 


CM 

k, B, C 



(1 

= 1 - 

■ 4) 


(1 

= 5 ■ 

- 20) 


A2, 





^8 



“2i 

£.3 


l|M|l 





0 

O' 

Q 


R. R 




£b' S 


"a' “a 


‘'^B' “b 

2 ^; ( 1 = 1 - 20 ) 


( 1 = 1 - 20 ) 


Center of mass of the satellite 
Satellite subassemblies 

Masses of rigid bodies having moments of Inertia 
Point masses (scalar) 

Orthogonal unit vectors, fixed In the body A 
Orthogonal unit vectors, fixed In the body B 
Inertlally fixed orthogonal unit vectors 
Total mass of the satellite (scalar) 

Nominal location of the center of mass of the body B 
Origin of the Inertlally fixed coordinate axes 
Reference point, fixed In the body B 
Position vector of Q relative to 0, and Its B-based 
matrix 

Position vector of center of mass of the body B and 
Its matrix, both based In the Inertlally fixed axes 
Displacement vector of the center of mass of the body 
B and Its B-based matrix 

Angular velocity vector of the body A and Its A-based 
matrix 

Angular velocity vector of the body B and Its B-based 
matrix 

Displacement vector of the masses m^ In B-based 
coordinates 

Displacement vector of the masses m^^ In A-based 
coordinates 
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Rotations of the rigid bodies In B-based coordinates 


^ 1.1 


’ ®1,2 ’ ®1,3 
(1=1-4) 


A®1,1 ’ A®i,2 * A®1,3 
(1=1-4) 


= Rotations of the rigid bodies In A-based coordinates 


1 . P 

I|m||^ and ||M||g 

(1 = 1 - 20 ) 

T, T 

T^. T (1 - 1-4) 

H 

(1 = 1-4) 

(1 = 1 - 20 ) 

(1 = 1-4) 


Force vector on satellite and its matrix, both based 
In the Inertlally fixed axes 

Total masses of the bodies A and B, respectively 

Force vector on and Its B-based matrix 

Torque vector on the satellite and Its B-based matrix 
Torque vector on the rigid bodies and Its B-based 
matrix 

Angular momentum vector of the satellite about CM 
Angular momentum vector of the rigid bodies about their 
centers of mass 
Nominal position of m^ 

Center of mass of m^ 

Position vector of relative to Q and Its B- 

based matrix 


N 


a 


1 


0 


= Acceleration of masses in the inertlally fixed 

axes 

= Matrix for transformation from Inertlally fixed axes 
to B-based coordinates 


N 


(■) 


= Inertial time derivative 
= Time derivative In body fixed coordinates 
= Skew symmetric matrix operator 

= Displacement vector of beams In B-based coordinates 


( ) 

( ) 

(1 = 1 - 8 ) 
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Xi (i “ 1-4) 

bPi (1 “ 1-8) 

pP, (i - 1-4) 

S*’a ’ 8*’b 

*^ 1 . 1 ’ 1,2 * 1,3 

(1 = 1 - 4 ) 

(1 = 1-4) 


Displacement vector of plates In B-based coordinates 

Displacement vector of shell B In B-based coordinates 

Mass per unit length of beams 

Mass per unit area of plates 

Mass per unit area of shells A and B 

Local orthogonal coordinate axes for beams, and fixed 
w.r.t. the B-based coordinates 

Matrix for transformation of p^-axes to the B-based 
coordinates 

Inertia force on the 1**^ beam element In B-based 
coordinates 


(1 -=■ 1-4) 

b**! 


B^ • ^ 


•>*^ 1,2 »“^ 1,3 


b\ 


b^-1,2 b^l,3 


Inertia force on the 1^^ beam element In local 
coordinates 

Elastic deformation vector of the l^'' beam In local 
coordinates 

Position vectors of the 1^^ beam element from the 
reference end In local, A-based and B-based- coordinates, 
respectively 

Position vector to the reference end of the 1^^ beam 
In B-based and local coordinates, respectively 



Thermal bending moments of the 1 


th 


beam 


Modulus of elasticity for the 1^^ beam 

Moments of inertia of the 1^^ beam about local 
coordinate axes 
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K and K,^ 
i,2 ^i.3 


b'^l 


b’^i 


“i.2 “i.3 

* * 

K and K 
^1.2 '■1.3 


r^i.Jk 


p“l 


pi pi pi 


° Thermal curvatures of the 1^^ beam about local 
coordinate axes 

<= Thermal bending constant for the 1^^ beam 

“ Characteristic time for heat transfer across the 1^^ 
beam 

= Attitude angles of the l^*' beam w.r.t. the Sun 

= Maximum values of the thermal curvatures of the 1^^ 
beam 

c Mass moment of Inertia matrix of the 1^^ rigid body In 
B-based coordinates 

= Stiffness of the 1^^ plate 

= Modulus of elasticity, thickness and Poisson's ratio 
of the 1*"^ plate 


pT^ , = Temperature distribution and thermal coefficient of 

expansion of the 1^^ plate 


k. 
P 1 




p’^i.O 


Thermal constants for the l*'^ plate 


p'i 

* 

p'l 


^B.r ’ ^B.e ' 


s^B.l * s®'b,2 ’ s*’b,3 


b*'b 


Attitude of the Sun from the plate nominal normal vector 

Flexural change of attitude of plate element from 
the nominal normal vector 

Radial, tangential and axial deformation of an 
element of shell B In B-based coordinates 
= Inertial force components per unit area of shell B 
in B-based coordinates 
Thickness of shell B 

Nominal radius of shell B and plates 1 and 2 
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8^B ’ s^B 


T 

8 B 


8®B 


* 

8®B 


= Poi88on'8 ratio and the modulua of elaatlclty of 
ahell B 

•> Temperature distribution of shell B 

» Attitude of the Sun from the nominal normal of an 
element of shell B 

- Flexural change of attitude of the normal of an 
element of shell B 


s'S * s'^B * s'^B.O 

Jl 

S' % 

Pi' S.l' jS.k 

u 

Xj» Aj 


J = 1-3 

p*^,j' J “ ^'2 
p'^.O 

fL , J = 1.2 
P 1 


Vi 

[T]^. 1 = 1-12 


Thermal constants for the shell B 

Identity matrix 
Defined by Eq. (1.10) 

Defined by. Eq. (1.3) 

a scalar, first zero of J o(x) 

scalars,, second zeros of Jj (x) and Jj(x), 

respectively 

Defined by Eq. (1.23) 

Thermal curvatures of the 1^^ plate 

k, . X, 

pi pi 

Thermal bending moments of the 1^^ plate 

Position vector of the 1^** plate element In B-based 
coordinates 

(3 X 3) square matrix of zeros except for 

[J = 1 

Laplaclan operator, defined by Eq. (1.27) 

(3 X 3) square matrices, defined by Eqs. (1.107) 
through (1.115) 
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ri,e 

s 

Plane polar coordinates fixed li 
’^l 

r 

B 


Xi.j 

S 

Defined by Eq. (1.34) 

J=l-15 

a 

Defined by Eq. (1.35) 

w(z) 

s 

Mapping, defined by Eq. (1.36) 

j = 4-6, 

a 

Matrices, defined by Eq. (1.39) 


I A ,], { q. } 
p i,k p’l 


r B, 

[Wrl 


f^l] 


s**3,0* s'%,z 


s'^B.O 


T 

8 B1 


b“b 


Matrix for tranaf omatlon from rectangular to 
cylindrical coordinates 

‘ f" £r-fi 

Position vector of an element of shell B In B-based 
coordinates 

(3 X 3) square matrix of zeros except for 

[ii]l,l=l 

Thermal curvature of the shell B in tangential 
and axial directions 


s'^B ■ s^B 


Shell thermal parameter, defined by Eq. (1.49) 
Thermal coefficient of expansion of shell B 


a,B 

b^, 1 = 1-11 
[XjLj], 1 = 7-9, 


° Angles, defined In Pig. 6 
= Coordinates, defined by Eq. (1.58) 

° Shell parameters, defined in Eq. (1.60) 
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. j = 11-14 

T^, 1 = 1.2 

T!, i = 1.2 


[X^]. 1 = 10-12. 

W‘ [/] 

'^l.l* ^ 

k^^2* ^ = 13-16 
[Xj^], 1 = 13-15. 

[X^], 1 = 16-18, 

[/i], 


Ag. 2 .J. 


'i'l, 'i>2> '('3 


01. 6^. 63 


AB 


F T 
■^C’ -IfiC 


= Rigid body coordinates, defined by Eq. (1.77) 

= Plate reaction forces on 

= Plate reaction forces on m 2 

= Diagonal elements of [ Ij . 1 ,] 

r 1 , JK 

= Inertia torque vector on m^^ 

= Resistive torque vector on m^^ 

= Point mass parameters, defined by Eq, (1.70) 

= Constants for the springs connecting m^ to plate No. 1 

= Constants for the dampers connecting m^ to plate No. 2 

= Spring-mass-damper system parameters, defined by 

Eq. (1.76) 

= Rigid body parameters, defined by Eq. (1.93) 

= Lengths of the rigid bodies in the ^2 
directions, respectively 

= Rotations of the body B relative to the Inertlally 
fixed coordinates 

° Rotations of the body A relative to the body B 

= Matrix for transforming vectors in the A-based 
coordinates to the B-based coordinates 

= Contact force and torque vectors exerted by the body 
C on the body B In B-based coordinates 
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Contact force and torque vectors exerted by the 
body C on the body A in A-based coordinates 
expressed In S-based coordinates 

IK 

II«IIb 

Ratio of the viscous and elastic modulli of the 
body C (Eq. (1.111)) 

The modulus of elasticity, moment of inertia of the 
cross-section and the length of the body C 
The modulus of viscosity of the body C 



AS^U>A/!£0 Co A 

ri.e)(/a£.£ 
3 l 3 



Figure 


£>/=■ 7“/v<r 


^^7-£^i-/r£r 


P£AT£- 3 




y se£f<^ < 

"Ko 


1 
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Figure 


MO. 
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Figure 



Figure 4 
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. Expressions for the Displacement of the Center of Mass 

The expressions for the displacement of the center of mass are required 
for obtaining the equations for each Individual structural element. But 
these expressions can only be obtained after the deformations of the Indi- 
vidual elements are known. However, it becomes easier to grasp the analysis 
that follows when these expressions In general form are available. 

The basic assumed configuration of the satellite Is shown In Figure 1. 
It consists of two composite bodies A and B, each containing point masses, 
rigid bodies, beams, plates and shells. These two bodies are connected by 
a flexible connector C. The shift of the CM Is shown by the vector £ In 
Figure 2. Let us consider the body B and find , the shift of Its center 
of mass In terms of the displacements of Its elements. The case for the 
body A can be obtained In an analogous manner. The mass of C is negligible. 
Using the B-based coordinates. It Is easily seen that 

-B ° [] I + n/ pPk + // sPb 

where 1 = 1,2, and 13-20; j = 1-4; k = 1,2. The Integrals are taken over 
the lengths or areas of the elements as applicable. The length of the beam 
Is E. To evaluate the beam Integrals from the beam displacements In local 
coordinates, we set 

Hi = [U®] •^,al (1.2) 

in Eq. (1.1). The local coordinate system for fhe beams is a right-handed 
system of Cartesian coordinates, with Its origin at the clamped end and 
^ axis along the length of the beam. For each beam the three components 
of are set up as follows: 
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b'Jl.l “ 



b'^1,2 

2*^1, 2^*^^'® 

s’*!, 2^*^^ > 

b‘*l,3 

iqi 3(t)*[exp(PiS)-l] + jqj^ 3(t)*s 

s’*!. 3^*^^ 

4 

1 

1 


■ f bPl’Ild**® 
0 

= bPj'tPj’ • / 



(1.3) 


f BiT 


■*■ ®2*2‘'j,2 ®3*3‘*j,2 

l^l‘l‘*j,3 ®2‘2‘*j,3 ■*■ ®3*3‘*j,3-J 


(1.4) 


Here obviously *= [ ^ (e^’^-1) - 2] ) 


Bj ■“ “Y Bj 


(1.5) 


The sum of the beam Integrals is obtained by evaluating [Vj] for each 
beam and then using Eqs. (1.4). 


The rigid bodies and the point masses In the spring-mass-damper systems 
are assumed to move only In the direction of b^ (see Fig. 1). The masses 
at the ends of the beams can move In all directions. The shell Is assumed 
to deform only radially. The plates are assumed to deform only laterally. 

So we can rewrite Eq. (1.1) as; 
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1 




ll«lt 


S7'yi7.2-^18'yi8,2-^19'yi9.2-^2o‘y20.2+ s^B ^B.r sine-a^dedz;) 

^7 yi7.3-%8 ^18. 3-^19 ^19.3-^20 y20.3-^iyi'^2y2-^13yi3'^u'yi4 
+ "15^15 ■" “16^16 ■" pPl // X/A + pP2 ;; X^dA) 


ll«lt 


I / bPj2d"= 


( 1 . 6 ) 


We have taken the shell displacements later on as 


^B.r “ bi(t)-4>i(e,ti) . 


(1.7) 


where <t'j^(9,zi) are purely spatial functions. Similarly, the plate displace- 
ments have been taken as 


Xi = “l,l^")’f'l,j^*'«®^ +1 j = 1 and 7-11. (1.8) 

15 * 

X 2 = “ 2 . 1 ^*^^ ■'•' 2 , 1 ^'^’®^ ■'■ ^ X 2 j(t)''i' 2 j('=^>®); j = 1 and 7-18 . (1.9) 


Equations (1.3) to (1.7) will become clear when we come to the individual 
analyses of the beams, plates and shells later. It Is evident from Eq. (1.6) 
that when the Integrations there are carried out, ^ becomes a linear 

“D 
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combination of only the time-dependent functions in (1.3), (1.7), (1.8), 

(1.9) and the ^(t), 1 = 1,2 and 13-20. 

* 

Now, let ^ column vector of all the time- 

dependent functions , 1,9^ • j » ®2 1 ’ ^1,1 ’ ^2 1 **1 “®“" 

* 

tinned so far. These generalized position coordinates define the 

deformations and displacements of the spacecraft. Then we can write ^ 

— B 

explicitly as 

Cj = [Gb1*% 

where [G„] is a (3 x HQ) matrix of known constants. The detailed forms 

O 

if 

of [Gg] and will be given in Appendix 1. 

2. Beams 

(a) Inertia Forces 

Referring to Figures 2 and 3, the equations for the inertia forces on a 
beam element are first obtained. We follow Likins [10] in doing this. 

Let dm^ be the mass of a beam element in body B. 

‘•“i = b^l * 

Let the position of the element with respect to the CM of body B be given 
by ^ + n. ) . Then the inertia force is given by 

D B 1. — 1 

b^^= + Cg + ii^ + 2(^ X (^ + n^) 

+ (0^ + (^ X X (^ + + ^8i + Hj)] 

where ^ locates the CM of the body E in inertlally fixed coordinates. 

Let V and oj be two vectors given in matrix form by 
T T 

V= [Vi, Vz, Vj] and a= [wi , (i)z,(i) 3 j . If we define 
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0 

-Vs 

V2 

V = 

Vs 

0 

-Vi 


_-Vs 

Vj 

0 

then V X CO e Vbo ® -SjV « *-co x v. 

Therefore 



+ C3 + + 2S^(C3 + 

+ (Sj + agSfjHCg + + 38^ + n^)] . (1.12) 


As mentioned befote, the axes of the local coordinate frame corresponding 

to the 1^** beam, are p 1 1 1^4 0 ^*4 •4* This la a right-handed system 

X}X Xy <9 

with U4 1 being always along the axis of the beam pointing away from the 
li-*- 

CM. p, , Is always parallel to the bs axis. The origin Is always 
l»d 

at the clamped end of the beam. 

So [p^] , the transformation matrix for converting B-based vectors into 
p-based vectors. Is given by 


r 

[Pl] 


'l,l’'’‘ 

^‘i,r*’^ '^1.1”” 



'l,2’‘» 


* 

(1.13) 





! the beam displacements in local 

coordinates, then 


b% 



(1.2a) 

°i 

= [ViI-CbSi) • 


(1.14) 


Then In local coordinates, the Eq. (1.12) becomes 

b^p^ = - “ d”ilvies + »'iS + b% + 

+ 2(p®-S3(p®)'^)3^ + (p®i3 + yXag)C3 + pJcSg + SSeSSeXv®)’' 

(1.15) 
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Let us define 


(y®)'£)g* (v®)^ = to* (1.16) 

and 

(U®)S'gSi' 3 (U®)’' = to** • (1.17) 


Then Eq. 



(1.15) becomes 

^p^-ds^.[y®0Xg + U®Ggq* + (lU^SjGj)^* 

■*■ 2“b ‘b^l H * ®i 


(1.18) 


+ b’l>l 


(b) ThermoelasClc Forces 

The thermal oscillations of the outstretched beams can cause considerable 
changes In the attitudes of the spacecraft. This problem was first success- 
fully modelled by Etkln [27J and later on by Yu [28]. Here we follow the 
method of Yu, and assume that the thermal curvature of the beam Is linearly 
proportional to the local heat input from the sun. 

Let K and K be the thermal curvatures about the y, , and 

^1,2 ^1,3 

y^^ 2 axes. Then we have the relation from Yu [28] as 




_ii2» - 
3t 


b'^1 


+ co8(a^ 2 + ) 


(1.19) 


for the y^ 2 direction and a similar relation for the y^^ 2 direction. 
Eq. (1.19) has the series solution (see Appendix 2) 


where 



b’^l 


02 


3 

3s^ b**!,!' 


bT^[cos a^^2 - - bV* + b^i® 


...)sln 
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If we set , k. •, T. « k_ * then for small values of , T. , we get 

DiDiT^2 

\ 2 = “1.2 ■ b<»i.2 - b^i • ^ “1.2I 

(1.20) 

and a similar equation about the )i .. axis. Then the thermal moment on the 

beam, , M_ Is given by 

‘>^.2 

The thermoelastic shear force gradient on the beam Is 


ii^bV = w<bV3> • • 


( 1 . 21 ) 


(c) Elastic Forces and the Equation of Motion 

If , f„ Is the elastic force gradient along the beam, then the equation 

D £ 


of motion Is given by 
3* 


“ ^.2 ^b“r^^ “ 3s j • 


3s‘ “ *1 


3s 

As 


The beams are assumed to be axially rigid, so 


b^E,l ° 


* .. ^b*^l,2^’ 


b E,2 'b r 'b 1,3 


3s’ 


and 


b^E,3 “ ^bV' Vl,2^ * ^b‘*l,3^ * 


so the final governing equations become 
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0 


^bV'^b’-l,2^['^^b‘*l,3^ '*' 3 ^b‘’i,3^|’ ®^““l,3 j 




• A **. 


•" \.l 


+ “bVS'IbI + “b •" “b > <! ^1.2 

^ 1.3 


0 . 


( 1 . 22 ) 


(d) Galerkln's Method 

The solution Is approached through Galerkln's method. The solution vec- 
tor Is assumed to be a finite sum of known space-dependent functions multiplied 
by unknown time-dependent functions. To keep the dimension of the reduced 
problem small, the solution Is taken to be an exponential flexure term 
superposed on rigid body rotation and translation terms. Thus 


b% 


^“>1,1 ' 




b‘^1,2 

> ° 

< 

iqi,2(t)*(e‘’‘^-l) + 211,2^*^^ + 3'*i,2('^^ » 

.b^i,3 . 



iqi,3(t)-(ePi«-l) + 2‘>1,3^*^^’® 3'll,3^‘^^ , 


^ . (1.3) 


When Eq. (1.3) Is substituted Into Eqs. (1.22), the left-hand side becomes 
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0 




The resulting equations become free of the partial differential operators. 

Each equation Is then multiplied by 1, and s, successively, and Integrated 

between (0,2). This gives us nine equations In seven unknowns. We reject the 
two equations generated from the first row equation by and 1. This 

leaves seven equations In seven unknowns for each beam. These equations also 

A 

Involve the other 103 unknowns In the vector Another 103 equations will 

be generated from the other structural elements. The vector 0Xg Is ex- 

* 

pressed later as a function of 

If we form a (7 x 1) vector for each beam, given by 




^‘’l,l’l‘*l,2* 2‘'l,2* 3‘'l,2’ l'*l,3’ 2^1,3’ 3‘’i,3^ 


then the generated equations for each beam will be of the form 


^b^ll^^b^l^ ^b*'12^^b‘*l^ ^b*^13^^b‘*l^ 


where ° constant square matrix. Also, Is a constant. All 

other quantities are functions of ^ and . The reduction of Eqs. (1.22) 
to Eqs. (1.23) In the cases of the four beams (Figure 1) Is shown In Appendix 

3. 

The gravity gradient and other environmental torques are not considered 
to be acting on the beams. They are Included In the equations for the 
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composite bodies as concentrated torques. The continuity conditions at the 
free ends will be used In the equations of motion of the tip-masses. The 
conditions for the shell-beam Junctions will be obtained later. 


3. Plates 

(a) Load System 

The load system on the plates Is shown In Figures 1 and 4. As In the case 
of the beams, the general equations of motion for the plates under Inertial and 
thermal loads will be derived first. These will be reduced to time-dependent, 
second-order, ordinary differential equations by using Galerkln's method. 

The continuity conditions at the Junctions to the spring-mass-damper systems 
and to the rigid bodies will be used In the equations of motion of those sys- 
tems. The conditions at the shell-plate connections will be defined later. 


(b) Inertia Forces 

To keep the governing equations linear In the space coordinates, the In- 
place deformations are considered to be negligible. The lateral deformation of 
the l^** plate, 1® parallel to the axis This makes the Introduction 

of a local coordinate frame unnecessary. 


Let be the matrix 


0 0 0 
0 0 0 
0 0 1 


Then the Inertia force on a 


plate element Is given by 


p^i “ - pPit'^alK + Xi) 

+ (“b + V + (<^ + , 

where Is the position coordinate of the element of area dA 

coordinates. Equation (1.23) Is only a scalar equation. 


(1.23) 
In B-based 
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(c) Thermoelastlc Forces 

The distribution of the thermal bending moments over the surface of the 
plate is assumed to be first-order, time-dependent. This model is obtained in 
a way similar to that used in the beam. The plates are assmned to be thin and 
homogeneous, so that a linear temperature distribution across the thickness 
of the plate can be assumed. This assumption also makes the thermal bending 
moments and curvatures at any point on the plate in two orthogonal directions 
equal. 


Let , and k, . be the thermal curvature of the plate parallel 

p 1,1 p 1,2 

to the plane formed by axes b, and b». Also, let i ■> 

*Ti. P^*^ 

be defined as the average curvature of the 1*'*' plate. Then is an in- 

variant of the plate with respect to a rotation of coordinates in the plane 
of the plate. Now let k , be given by 


k cos( B, + B, j) 

pi pi p 1>1 


where 




the angle between the plate nominal normal and the sun vector 
the rotation of the plate surface normal, due to flexure, along 


The solution of this equation (see Appendix 2) is given by 


p 1 p i pi p 1,1 pi P 1,1 pi P 1,1 


3 '«* 

■ p""^! ' p®i,l ■*’ 


Keeping only the first power of 


T, , the solution takes the form 
P i 


131 



similarly we have the solution along ^2 

= p'^i.0l'^°« p\ - (/t,2 - • p\.2)sl” pPil 


(1.25) 


where 


n = k . • T. = the maximum value of »c. , 

p 1,0 p 1 p 1 p 1,1 

(C, , = a constant. 

P 1>2 


Let xj^ and X 2 be the distances In the directions bj^ and ^2 • Then 


* ^Xl * 

p^l,l " 3^ • p®i,2 “ ^ 




and 


3Xi 3Xi 

r-. - ( K. , + K. ,) » (2 K. -cos e. - K. „[(-5-i+ TS-i) 
pi p 1,1 p 1,2' p 1,0 p'^l p 1,0‘ 3xi 8 x 2 ' 


3 3Xj 

- T. + ■5^)]sin ej 

p 1 3t'3xi 3x2 pi 


(1.26) 


If ^L and are the bending moments along the b. and b, 

P 1,1 P ‘■1,2 

axes, respectively, then 

P\,l ° P^i ' P’'i.2^ 

and 

M_ = D, [ K, , + p, • I ] • 

P\^2 P 1 P 1.2 pi P 1.1 

It Is assumed that the thermal twisting moments In the plate are absent. 
Thus the equation of motion Is obtained as 


<A.i * <p».,2 * p\,,> - 


) = - q 


where q is the plate loading, M, , , M, . and M, ,, are the usual 

P 1.1 P 1.2 P 1,12 

elastic bending and twisting moments of the plate. The above equation can be 
written as 
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- p®l • + p^'l>’”p'^l.l + p'^i.2^ = - ■> > <1-27) 

where 

j 

3x| 3x| 

a^d pD^ and are the stiffness and Poisson's ratio of the plate, re- 

spectively. So from (1.23), (1.26) and (1.27), we obtain 


V? = + 


+ (1 + pVp-'l.O 




+ pPi[^3][0Xg + ^ ^ 


+ («b + %VVb-*-4-*-pV “ °- 


(1.28) 


In Eq. (1.28), Qx^ Is the Inertial acceleration of the body B In B-based 
coordinates. As was done in the case of beams, SX^ Is replaced by 


vi^iib 


=n^^Bc 


(1.29) 


where Is the contact force between the bodies B and C. Is ex- 

pressed In turn as 

2fic = 4 ^1" ^2 >Vb -^ 3 ^ + ■^7 - - ^9^ • (^-29) 

Equation (1.29) will be explained In detail later. 


(d) Galerkin's Functions for Plate Nos. 1 and 4 

Before applying Galerkin's method to Eq. (1.28), It Is first transformed 
into plane polar coordinates. Let 
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XI = ri COS0 and X2 = ri alnS 


where 0 < ri < a . 

— — o 


Now, with 


v! = 


_3f_ ^ JL_ _3_ 
3rf n 3ri 


ri 30^ 


3 

3x1 


COS0 


3 

3ri 


sln6 


3 


3 

3x2 


3 1 

COS0 -s sln6 

3ri n 


3 


we have 


3 

3x1 



(sln0 + CO80) + — ~ ) 

3ri n 3 Q 


So the left-hand side of Eq. (1.28) Is transformed Into 


+ {(1 + yj)' <i_^o peiHsins + CO80)(^ + ^ 

do 

- + pVVi.o ■ p^i pSiXsine + • 

Now let a^r = ri, so that 0 ^ t ^ 1, and 0 ^ 0 ^ 2 ti. Then V{ = ■ 


B 


z 

where V 


Expression (1.30) becomes 


3r 


2 ^ r 3r 


1_ _3f_ 

r^ 30^ 


- {—(1 + y,)- K, . 
a' pip 1.0 


p^'l^^l.O pei}(8in0 + co80)(^ + 

• p^i pPiX^i^e + e°«0) + 7 • 


(1.30) 
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Also 


Since Is a vector [0, 0, Xj]^, = 0 = • 

“ “ ^“b 1 ■'■ “b 2^^ ■ vector 

T 

[a^r cos6, a_r sinS, R . ] . Then 

O . D 1- 

-«3(Sg + SgVA “ ^“b,i“b.3 ' %,2>V 

%,2“b, 3 ^ ^.l^V - K.l + “l.2>^ • 

So the right-hand side of Eq. (1.28) becomes 


-pPll)(l - + <‘^,i“b.3 - ^.2>V 

+ <“b,2“b, 3 + '^.l^V - Kl ^ 

+ pP^[-^ 3][Ggq* + {2S!g + ^ (i - T2)}Ggq* 

+ {Sg + + [^ (i - Tg)}G^q* - f] 4 r^ (Tgi + T 36 ) ] ( 1 . 

Equating the expressions (1.31) and (1.32), the equation for the plates is 
obtained as 

p°i p; + 1^(1 + pUi)p<i,o p6ij(^i"e + 

L“b IfB J 

- p^l + p^l^^l.o p6i|(si“e + cos6)(^ + i ^)V^Xi ] 

^ pPi'Xi - K.l - ^.2> V 


^^.2“b, 3 + ‘^,1^V - K.l 

+ pPi['^ 3 ][Ggq* + {223 + Tj - + 
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0 . 




+ T^e )] = 


Now the choice of the Galerkln's functions must be made so that the solu- 
tion matches the load system and satisfies the boundary conditions. As there 
are four spring-mass-damper systems and a central junction on the plates 
(numbers 1 and 4) , the solution should contain five terms which are singular 
at the connecting points. This has to be done because the loads are taken to 
be concentrated forces. We also apply the boundary condition that the space- 
dependent solutions be zero on the boundary r ° 1 of the plates. It Is also 
asstmied that the nature of the solution Is close to the free vibration mode 
for thin plates. Therefore the solution Is taken to be 


^1 ° + Xi^4(r,t)slne 

+ Xj^^ 5 (r,t)cos 26 + g(r,t)sln 26 + r*log r 

+ Xi.u(t)^,^ll(r.6) . (1.34) 

The functions xji (r,0) are basically the static point load solutions for 
1» J 

plates. 

Xl 2’ Xi^ 3 > Xi 4 > Xi 5 Xj 8 chosen In the following way: 

Xi 5(r.t) = ^ 

“ \,4'^2^X2r) + 5(r-l)r^ + a^^g(r-l)rS 

Xi sCr.t) = a^_ 7 Jj^(X^r) + a^^ g(r-l)r + a^ ^(r-Dr* ^ (1.35) 

Xi^4(r.t) = 10-^1 (Xj^r) + a^^j^j^(r-l)r + 
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= ai,i3Jo(V^ + ' 

Here Jo(x). Ji(x) and J2W are the ordinary Bessel functions. 

Ao, Ai and A2 are the first zero of Jo(x), and the second zeros of Ji(x) 
and JzW, respectively. Apart from a quadratic term in r, the coefficient 
of 7^*"^ solution to the problem of a plate with a steady concen- 

trated force at the center. 


(e) Coefficient of g(t) 

The functional form of .(r,6), which have coefficients X. 4» 

J = 8-11, are similar. The expression for \l>. „ is derived here. The 

ifO 

forms of 9» 30 ’I'j obtained in a similar manner. 

Let g be the mode corresponding to the spring-mass-damper system 
containing m^^^. Let z = r exp (16) represent the points in a complex 
plane. Let cXg = (Og + iOg be the coordinate of the load point on the 

plate number 1, corresponding to mj^g. The plate geometry is the unit circle 
r = 1. Let the function w w(z) map this unit circle onto another unit 
circle w = 1 and the point z =■ a to w = 0. This mapping is given by 

z - Oo 

w = (1.36) 


1 - 

when Og ■» (Og ^ - iOg “ bbe complex conjugate of a. Then the solution 

to the problem of a circular plate with a concentrated load at z = a in the 

z-plane corresponds to the problem of a circular plate with a concentrated 

load at its center in the w-plane. The problem then reduces to choosing a 

blharmonlc function i|). „(w) which 1) will be zero at w = 1, li) will have 
1,B 

a logarithmic singularity at w ° 0, ill) will not be harmonic, Iv) will tend 
to r* log r as Og 0, and v) will be real valued. It can be noted here 


that 
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and 


V* S 4 


V' = 16 


3z3z 


3z^3z' 


This leads us to choose the required function from the list of functions 

■k 

given below. 


* 

>('1 


* 


* 

lj)3 


* 


<11* 


where 


(z-a ) (z-a )log(ww) 
o o 


(z-UgJw log w + (z-a g)w log w 


Re[(z-a )(z-a )log w] 
8 8 


Re[(z-a ) (z-a )log w] 
8 8 


Re((z-ag)w log w] 


Re((z-ag)w log w] 


z - ag 


1 - agz 


It can be shown (Appendix 4) that <{ii Is the most suitable choice. Thus 

g(r,0) ■= (z-Og) (z-oig)log(ww) • (1.37) 

Let Im[a„] = 0. Then a = a = a real number. Therefore 
o o o 


o o 

/ r^-2a r cos8 + a^ 

ip = (r^- 2a r cos9 + a*)log ( ^ 

® ® n-2a„r cose + air* 


(1.38) 


(f) Final Equations for Plates 1 and 4 

Equations (1.34), (1.35) and (1.38) are substituted into Eq. (1.33). 

The resulting equation is then successively multiplied by the 21 spatial 

functions introduced In Eqs. (1.34), (1.35) and (1.38), and then Integrated 
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between the limits of 0^6^ 2rr, 0^ r ^ 1. This leads to a set of 21 
second-order ordinary differential equations of the form 


+ ^SS^B + = ^pV 


(1.39) 


r ’*1 

In this equation, ( q./ la a (21 x i) vector consisting of the unknown 
P 1 

time-dependent plate functions. The detailed form of Eq, (1.39) Is shown 
In Appendix S. 

The plates 1 and 4 have similar equations with the same number of un- 

r *•» 

known time-dependent functions • 


(g) Equations for Plates 2 and 3 

The basic equation for the plates 2 and 3 are also given by Eq. (1.33), 

but the load conditions are different because there are two rigid bodies 

attached to each of these plates along with the spring-mass-damper systems. 

While the central load Is absent In these plates, each rigid body is assumed 

to be fixed at four points of the plates, so four concentrated load functions 

for each rigid body are added to the assumed plate displacements. The form 

of these functions Is similar to that of ’p. o- Therefore the solution to 

1,0 

the plates 2 and 3 Is taken to be 


^1 “ ^1,1^*^^ Xj 3(r,t)cose + 4(r.t)sine 


12 

+ Xi^5(r,t)cos 26 + Xi^e(r.t)sln 20 + Xi^g+j ^'^>’^l,6+j • 


(1.40) 


The functions X, ,, j=Z-6, are given by Eq. (1.35), and the functions 

1 ,-j 

j, j = 7 - 18, are given by Eq. (1.38) where Og Is replaced by the 
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corresponding coordinates of the load points. 

The final equations for the plates 2 and 3 are also similar to Eq. 

* 

(1.39), but In these cases the generalized coordinates ® ^ 

vector. The detailed form of these equations Is given In Appendix 6. 


4. Shells 

The analysis for the shells A and B are similar. Only the procedure 
for the shell B Is given here. 

The shell B Is assumed to be a uniform, thin. Isotropic, circular 
cylindrical shell. For Che elastic analysis, the linear equations of Vlasov 
[36] are used. The analysis of the thermal effects follows that made by 
Kraus [37]. The orientation of the coordinate axes Is shown in Figure 5. 


(a) Inertia Forces 

T 

Let ,, F_ ,] be the Inertia force vector on a 

& li 8 uyX S D f ^ 8 D}.} 

shell element In the B-based rectangular coordinates. 

Let F^ = [ F„ , F„ F„ ] be the Inertia force vector on the 

sHB s B,r’ s B,6' s B,z' 

shell element of area dA, In the B-based cylindrical coordinates. 
Proceeding as In the case of beams and plates, we have 


S^B= S ^B + V + (“b ^ + ^B>J 

(1.41) 


where 


Rjj = position vector of the undeformed shell element In B-based 
coordinates , 

T 

1 . r.] = displacement vector of the shell element 

D,/ Up J 

In rectangular B-based coordinates, and 

d = p • dA . 
m s B 
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r T 

Let 5- “ [C_ > ?„ a. ] be the displacement vector of the shell ele- 

-‘fi B,r D,y B,z ^ 

ment In B-based cylindrical coordinates, and 




CO80 sln6 0 

-sln0 COS0 0 

0 0 1 


Then 


and 


Let 


and 


Then 










{^F^}= - dm[u^0^ + + % + 2w'i2^C3 + 2u^ES^(y^) 5^ 

-•*T *T •* T 

= ^ <“B.r ^ %,rA 

In this analysis, It Is assumed that 


So 


“b.t^b 


^B,0 “ ^B.z “ ° 


B,3 


(u_ ,8ln0 - i0„ _CO80 
B,1 JS,Z 


■’B.r 


(1.42) 
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and 


** I 


-((ijg + “g j^cos0)((Ug jSin$ - oj^ 

1 “b.3%,2«^“® + “b.i“«®) 


Kr ■ 


With 


we have 


{Rg} = [agCoaS, agSine, z]' 


ri 0 0 ■ 
0 0 0 
Lo 0 0 


where 


l^l] 


Let {0^ '''S''’ B ''' ' Then, using the 


results obtained above, F„ is given by 
* s B,r ® 


sS.r “ " K.3 ■'■ -S,2“®®>"JS,r 

-'■ z(Uj^2‘^°°®““b 1®^"®^ “ ®B^S "S 2*^°®®^* ■ 3 

+ cosO'fi + 8 ln 8 ’f 2 } . (1.43) 


(b) Thermoelastlc Forces 

The distribution of the thermal bending moments over the surface of 
the shell is obtained in the way used for plates. The shell is assumed to 
be thin, homogeneous and isotropic, so a linear temperature distribution 
across the thickness of the shell can be assumed. This assumption also 
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makes the thermal bending moments and curvatures at any point on the shell 
In two orthogonal directions equal. 

Let o and K_ be the thermal curvatures of the shell In the 

sB^v SByZ 

tangential and axial directions. Then 


■ ■ . 8 + a , B,e ^ k_( 6 + B ) 

3t sS's’^B ^ s'^B.e^ 


(1.44) 


where 


= the angle between the shell nominal normal and the sun vector at 
any point, 

* 

A the rotation of the shell surface normal due to flexure In the tan- 
s 

gentlal direction. 

As shown before, keeping only the first power of , the solution to 

Bq. (1.44) takes the form 

s'^B,0 “ s’'b, 0^®°® 8®B " ^s®B,e ■ s'^B * s^B,6^®^“ s®B^ • d-^S) 

Similarly 

K „[cos 6_ - ( B* _ - T * B* )sln B„] (1.46) 

6 ByZ 8 ByU SB 6 DfZ 8 B 8 B^Z 8 B 


where 

b!! “ the rotation of the shell surface normal due to flexure In the 

s B,z 

axial direction. 


Now from Eq. (2.27) of Kraus [37], 


2- s“t„ = s®B - s''b><sS, 6-^ sS,z> 


s“t “ 24 ‘ 8®b'^s*'b^ ^8^,6 b''b,z^ 

B 


(1.47 


Let C ° + defined as the cylindrical surfaces of the shell. It 

Is assumed that the temperature distribution across the thickness of the 
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shell be given by T (6,z,?). Then, combining Eqs. (2.28) and (5.6) of 
Kraus [37], we have 


= s^B • s“b . 

Let be defined such that 

“ 12 ‘a®B ■ s“b ' s''b ‘ s''^Bl * 
From Eqs. (1.47) and (1.48), we obtain 


(1.48) 


( o + <» .) • 


(1.49) 


a B1 2(^0^) '■a B,6 ^ a B,z 
This result will be used to derive the equation of motion of the shell. 


(c) Equation of Motion 
Let 


“b^i 


where z Is the distance along the ^ -axis. We also define 


= 


3zf 


36* 


The equation for the radial displacement of the shell, obtained by com- 
bining Eq. (13.2), Part 2, of Vlasov [36] and Eqs. (6. lid) and (6.13c) of 
Kraus [37], is given by 


3" 


h * h* 

^ (V- + 2V* + DV"? - ^(1- u )( ^ - 
12a* 6a* ® ® 3z![ 3z!38“ 

D O 


■)V*5,. 


^ - s>^B>^ ^B.r 


aS . s«>B " « 




where q Is the shell radial loading. In this analysis, q consists of 
only Inertia loading given by 
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(1.51) 


- = s^B.r • 

Therefore from Eqs. (1.51) and (1.43), we have 


“ „P, 


’b|^ ^B,r ■ ^“b,3 2^*^,! “b,2^^^ ^B,i 


^I'K.l - 26 + 26][V'53^^ - BV^Cj,,, 


«B.r>l - "t<“B.l - ■4.2>«^" 26 - 2co3^^o>3 ^cos 26] 


'tW '^"^B.r ■ ^B.r^ ■*■ 


Now, from Figure 6, It Is seen that 


cos ^6g “ slnB cos (big jt + 6 + a) (1.53) 

where a is the angle between the plane containing and and the 

plane containing ^3 and the sun vector, at t ° 0. For simplicity, a 
Is taken to be zero. 

Then, from Eq. (1.46) and Eq. (1.53), we have 


^‘(sS,z> " - sVo 


“ s*'b,0^*^sS,z " s^B * 8 ®B,z^®^“ 8®B^ 


- “-<“b, 3‘ + 9) + s®B^‘^s®B,z " s^B + s®B,z>l- 


Similarly 


^‘^s'S, 6 ) - s'S,o£®^“® “«(“b, 3 *^ + e) + sin 3 eBV'(^B *_0 - ^Tg • ^Bj_g)]. 
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Hence, from Eq. (1.49), we obtain 

“«(“b, 3‘ + ®> + I sin ^6g(l - -^) 


•f 


= 8lne[cos(aij t + 6) + -^ s1“(“b, 3'^ + 6) d - Jg’ 

SB * B * 




(1.54) 


Now substituting Eqs. (1.52) and (1.54) In Eq. (1.50), we obtain the equa- 
tion of motion of the shell as 


8 5 /n** . 


12a: 


(V- + 2V^ + 1)V“£ - 5_|(i u 


, a>* 

+ r + -S: 

3zJ 8^ 


6a^ “ " 3z! 3zl38 

>B \ 




1“b,3 + lK,l ^ |lK,l - ‘^,2>-« 26 


^ 2“b,i“b, 2^^“ 26][V'e^^^ - - 16(^+ 5,^^] 

- ^I(“b,1 - "b,2>«^" 28 - 2<s^^^a,^^^cos - 4^3^^] 


+ f 1COS0 + fzslnsL + -^ (1 + y ) h|8ln6[cos(co ,t + 6) 
r 11 8 o B D B^J 
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(d) Galer kin's Functions 

The Galerkln's functions for the shell are chosen so that the solution 
remains compatible with the load conditions. The shell loading consists of 

1) the moments applied at the plate-shell Junctions, 

11) concentrated forces and moments applied at the beam-shell junctions, 

and 

111) the Inertia and the thermal forces. 

To accommodate the loading 1) , the solution must Include some exponen- 
tial functions of z i . To accommodate the loading 11) , the solution must 
contain the delta functions 6(zi-p^, 6 - 6^), where zi = and 
6=0^ are the coordinates of the junction of the 1^^ beam and the shell. 
The forces In 111) are taken Into consideration by the equation of motion 
Itself. The flexural displacements must be equal to zero at the boundaries 
zi = -p' and zi = +p" of the shell. 

Now, two expressions for 6(x-Q, given by Eqs. (4.7), (4,43) of 
Stakgold [38] are 

CO 

1 r 2 mtx mi5 

5(x-C) = T + I I ~ ~T' 

n=l 

0= (1.56) 

r 2 . iiTTx . mr£ 

“ \ 2 I 

n=l 

where £ Is the domain of the one-dimensional delta function. 

Since 

6[Z1. 6; p^, 0^] = 6(zi-p^)6(0-6^) (1.57) 

the form of the solution of the shell problem should be given by 
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^ = bi(t) + (zi+p')(zi-p")[b2(t)eP^^‘ + b3(t)eP^^‘ + ...] 

b,r 


. V V r, 1 . mTT(p'+pi) . mn(p'+zi) n n, 

+ ) ) [b., (t;sln — — ' sin cos noi cos nol 

m=l n=l U h 


+ I I [^ 3 +ij p(b) sin sin sin p 02 sin p6] 


m=l p=l 


i.1 




+ I I „(t) Sin — . 

m=l q=l *■* 


mTT(p'+z3) . uni(p'+zi) „ 

sin 7 - COS q0 3 cos q0] 


Ij 


V c r. 4 ^.> . nnr(p'+p 4 ) . mTi(p'+zi) . « . ni 

tbj^ (t) sin sin si„ ,64 sin r0] 

m®l r— 1 

(1.58) 

where li is the length of the shell. The last four infinite sums corre- 
spond to the four shell-beam junctions. Here 0i =0, 02 = Tf/2, 03 = tt 

and 04 = 3 y . Thus the cosine series is used for 6(6-0i) and 6(0-02). 
The right-hand side of Eq. (1.58) is then truncated to keep only the first 
and the second harmonic terms, and so the assumed solution is finally given 
by 

C- ^ = bi(t) + (zi+p')(zi-p")[b 2 (t)eP^'=* + b3(t)eP’*'*] 
b,r 

+ b 4 (t)sln COS0 + bs(t) sin cos0 


+ be(t) sin sin0 + b?(t) sin ~*-zi) 


+ bs(t) sin cos 20 + bs(t) sin cos 20 


+ bio(t) sin sin 20 + bii(t) sin sin 20 . (1.59) 

The coefficients p 2 and p 3 correspond to the free vihratlon modes of 
the shell. 
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Equation (1.59) is now substituted in Eq. (1,55) and the Galerkin’s 
method is applied. This gives eleven time-dependent equations of the form: 

+ [^7]{qg} + [XsHqg) + [X9]{qg} = {gAi,} * (1.60) 

- *, 

The vector is the (11 x i) vector with b^, i = 1 - 11, as Its 

elements. The details of Eq. (1.60) are in Appendix 6 . 

5. Beam-end Masses 

(a) Inertia forces 

The expressions for the beam-end masses of the body B are derived here. 

Similar expressions can be obtained when the body A is considered. These 

masses are assimed to be point objects. 

Let 1. = [Fj 1 . Fj oi F ,]^ be the Inertia forces on the i^*' mass 
i,i ^ f ^ i>"i 

along the beam local coordinate axes. 

Let ^ be the displacement vector of the i*"^ mass in B-based coordi- 
nates. Then, as before, 

F^ = -<n,[U®][0^ + Cg + 71 + 2£0g(C34i) 

+ (Ug + SgUg) (Cg + y^^ + Rj + r^)] . (1.61) 

Here (R^+r^) is the position vector of in the B-based coordinates 

in the undeformed state. Also, as noted before, 0 Xg and Cg are expressed 

r *1 

in terms of the generalized position vector I q_l . 

D 

(b) Elastic Forces 

Let = [F^ 2 » 3 !^ I’® elastic forces on the 1*"^ mass 

along the beam local coordinate axes. Then 
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F.' -1 = (Radial force exerted by the shell on the beam at the corre- 
spondlng junction point) + (the total inertia force of the 
beam) (1.62) 

2 ~ T'^®“sverse shear force along exerted by the 1^^ beam on 

the tip-mass 

2 ” Transverse shear force along exerted by the 1*'^ beam on 

the tip-mass. 

An expression for the shell outward load q is obtained from Eq. (6.11c), 


Kraus [37]. The expression is 


h* 

— V"? +5 

2 B.r B,r 




h^ 
s B 


12a 


B 


^‘•i-^(1 + sVsVVb.1 


s B s B 

where is the outward load on the shell at the junction with the 1 

beam. The above equation is rewritten as 


th 


E • h E • h E • h a 

s^B s“b . SB SB , s B SB . „2, „ V 

12(l-^PB)^a; 'B.t ^,.^^ 2^,2 ^B.r <s Bl> 


Zl = Pj 
■0=0, 


(1.63) 

Now the expressions for and T„, are substituted in Eq. (1.63) 

OyT S 15X 

and evaluated at give q^ in the form 

11 12 


‘'l = '=ij'’j ^ jJl2 %j^j • 

All c. . are constants. The details are In Appendix 7. 
l» j 

th * 

The total Inertia force of the 1 beam, , Is given by 

0 

* At r * • •* ** 

b^i = -['^'1 / (bPi)tb'li + (“b + “b 


(1.64) 


+ [pB] [0Xg + + (i^ + “b^)Cb]Msi 


(1.65) 
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where 


Ml] = 

1 

0 

0 



0 

0 

0 

• 



0 

0 

0 



The transverse shear forces along 

ik 

and are given by 


^1,2 “ ^b^l^^b’'l,3 

Pi 

£?!*• 

l‘*l,2^ 

(1.66) 

®’i,3 ° ^bV^b^i,2 

pi 

£?!*■ 

l‘*l,3^ • 

(1.67) 


Then ^ Is obtained from Eqs. (1.62), (1.64), (1.65), (1.66) and (1.67). 


(c) Equations of Motion 

The equations of motion of the beam-end masses are obtained by setting 


F. + F! = 0 . 

— 1 — 1 


(1.68) 


Rewriting Eq. (1.68), we have 


mi[Ui] [OXg + G^qg + y^ + 2o>g(Ggqg + y^ 

+ + yi + «1 + ■'i^i 


"jIi'i.iN * * 'i“» -'""‘’I 


^b®i^^b^l,3 l**!,!^ 

■ ^bV^b^i,3 P* ®*’* l**i,3 


(1.69) 
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There are three equations In Eq. (1.69) for the three variables y^ 2 

and y The four beam-end masses will be described by 12 equations 

1, J 

which will be of the form 

+ [Xlo]{q*} + [Xij]{q*} + [Xi 2 ]{qg} = {^A^}. (1.70) 

r *1 

Here ly ) Is the (12 x 1) vector of the displacements of the four beam- 
end masses. The details are given In Appendix 8. 


6. Spring-mass-damper Systems 

The masses In the spring-mass-damper systems are assumed to move only 
along the ^ -axis, so the displacement vector In B-based coordinates Is 
given by 


2i ” [0. 0. y^] 


Let 


[-»3] 


0 0 0 
0 0 0 
0 0 1 


(1.71) 


(1.72) 


Then the Inertia force on the masses Is given by 


= -m^[^3] I + 2<2b ^ + (“b + + S’b^ 

ot 




(1.73) 


Let ^ be the stiffness of the spring connecting the mass m^^ to the 
plate No. 1, and 2 X>e the coefficient of the damper connecting m^^ to 
the plate no. 2. Then the equation of motion of m^^ Is given by 

^1 “ *^1,1^^! " ‘‘l,2^yi ■ ^2^“l^^ (1-74) 
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where Is the position coordinates. In the plane of the plates, of the 

mass m^. Substituting Eqs. (1.71), (1.22), and (1.73) In (1.74), we obtain 
the equation of motion of m^ as 


“if^l - ^ “B,2>yi ^ %,1 “b,2“b,3>«1,2 <“b,i“b,3 “ “b,2>\,1 

■ K,1 + “b.2^'^1,3 ^“b ^ “b 

+ [-OajeXj] + - Xi(“i)] + ki,2^yi ' X2(V^ “ ° * 

Equation (1.75) Is a scalar equation. For each of the bodies A and B, 
a system of four such equations are obtained, corresponding to the four 
spring-mass-damper systems In each body. Let 

{ d<l } “ tyi3.yi4.yi5.yi5l 

be the (4 x 1 ) vector representing the motion of the sprlng-mass-damper 
systems. Then the four equations of each body is obtained In the form 

td^l^tj^ ^ td*2^tjq 1 + 1 


t^sltlB t\4l*lB t^lsl'^B 


taA.l 

d 4 


Details of Eq. (1.76) are given In Appendix 9. 


(1.76) 


7. Rigid Bodies 

The force and moment equations for the rigid bodies are now derived. 

As before, the required equations are derived for the rigid bodies In the 
body B only. The corresponding derivations for the body A are similar. The 
positions of rigid bodies are shown in Fig. 7. 
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(a) Force Equations 

The rigid bodies are assumed to be so constrained as to have a transla- 
tion only along the have two rotations about the ^ 2 “ ^2~ 

axes. So the motion of the rigid bodies with masses m^ and m 2 are 
described by 

= [y^. ^ (i.77) 

The mass m^ Is fixed to the plate no. 2 at the four points with position 
coordinates 1 = 7 - 10, In the plane of the plate. Similarly, a^, 

1 = 11 - 14 are the Junction points of m 2 with the plate no. 2. Let the 
reaction forces of the plate on mj^ at the points a^, 1 = 7- 10, be given 
by Ppi, and respectively. The reaction forces on 

the mass m 2 are given as ^2^2^’ ^2^2^ ^2^2^ ’ 

Summing the Inertia forces and the reaction forces to zero the force 
equation of motion for m^ Is obtained as 


mfCyi - - (‘^,2 ‘ “b,i‘^,3^''i,1 * “b,2“b,3>''i,2 


- K,1 “b,2>*'i.3 ^ + (Sg + 




( 8 ) 


,( 9 ) 


,( 10 ), 


(1.78) 


The plate reactions are obtained from the singular components of the 
assumed plate solution. The total shear force on a small circle of radius 
e at a distance r^ from the center of the plate due to the function 

fj = X 2 j(2 - oij)(z - aj)log(ww) (1.37) 


on the mass m^. Is given by 
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where 


Equation (1.79) Is derived In detail In Appendix 4. 
Substitution of Eq. (1.79) In Eq. (1.78) gives 


"l^yi -K.1'%.2 ^''i ■ (“b.2-“b.i“b.3^‘'i.1 ^^.1'^B.2“b.3^"i.2 

- ^ + 2iSg ^ 

10 

+ I'=<'3]0Xb> - 24 pDjJre^ f rj^(l-rj)X2 = 0 . 

The equation for Is similarly obtained as 

”2^y2 ■ ^“b,1'*^,2^^2 " ^“b,2"“b,i“b.3^®2,1 ^'^,1’*^,2“b,3^®2,2 

- K.1‘^.2^''2.3 + ^ + 2(3g ^ + (Sg + 

14 

+ [^ 310 X 3 } - 24 pD 2 ire^[ rj (l-rj)X 2 ^j] = 0 • 


In these equations, 3 ] Is given by 
[ J?3] ■= 


0 0 0 
0 0 0 
0 0 1 


(b) The Moment Equations 

Let Be the Inertia dyadic of the rigid body about 

troldal axes, and T. be the Inertia torque. Then 


where 


h “ - 

{§!> = ®i,2’ 


(1.79) 

(1.80) 


(1.81) 

(1.82) 

Its cen- 

(1.83) 

(1.84) 


155 



Equation (1.83) Is revritten as 


- ^ + 6^)1 - (“b + ®l) 

- r^i<^ - “b " r^l ' -i^ " r^l * ®i “ ®i " r^i • 




It Is assumed that form 


fr^ll 


r^l.l ° 


r^l,2 


r %3 J 


( 1 . 86 ) 


Then Eq. (1.85) reduces to 




^r^l,l%,l ■*■ ®1,1^ ■•■ “b,3^“b, 2 ■'■ ®l,2^^r^i,3 " r^l,2^^ 


" ■{ ^r^l,2^'S,2 ■*■ ®i,2^ ■*■ “b.3%,1 * ^i,P^r^l,l ~ rh.3^^ 

^r^i,3*^,3 ■*■ ^r^i,2 " r^i.l^ H, i“b,2 ^ “b,2®1,1 “b,1®1,2^^ 


(1.87) 


Let T' ^ and T' „ be the reaction torques exerted by the plate on 

i,X 

the rigid body about Che centroidal axes parallel to and 


respectively , 

. Then, from Figure 7, we have 


T * ts 

1.1 

1 . ( 9 ) ^ ( 10 ) 

^2.1 

- ■’a!! - 

(1.88) 

and 




t ' = 

1,2 

l(p(8) + p ( 9 ) 
2 '''^ 2,1 ^ ^ 2,1 

- • 

(1.89) 
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Slmi laxly. 


, ^ ifp(ll) + p(12) _ p(13) p(14).o 

^ 2,1 2 *''^ 2,2 ^ *^ 2,2 *^ 2,2 *^ 2.2 •'*'6 


and 


, ^ ifp(ll) + p(l«) p(12) _ p<13>so • 

^2,2 2'*^2,2 ^ *^2,2 *^2.2 2,2 ■' r 


(1.90) 


(1.91) 


Then the moment equations of the rigid bodies nos. 1 and 2 are given by 

t’^l.l* ^1.2* '^2,1> ^2.2l^ + in.l* n.2> ^2,r "2.2^^ = ° * 

Let {^q*} “ [?!• ®1 1> ®1 2’ ^2’ ®2 !• ®2 2^^ be the vector which 
defines the motion of the two rigid bodies. Then, combining Eqs. (1.81), 
(1.82) and (1,92), an equation of the form 

[^Ai]{r<i*} + ■*■ ^rS^^r***^ 


+ + [^7^93} + I^lgHqgJ ” ^r*4^ 

Is obtained. Details of this equation are shown In Appendix 10. 


(1.93) 


8 . Transformation Matrices 

The angular positions of the body B defined by three successive rota- 
tions of coordinate axes are shown in Figure 9. 

^2 S 3 Inertlally fixed coordinate axes. The B-based 

coordinate axes ^2 ^3 coincide originally with 2 ^^, X 2 and X^, 

respectively. First, a rotation by an angle about Xj^ brings the axes 

^ In the positions shown by D^. Second, a rotation by an angle about 

^2 takes the axes ^ Into the positions shown by E^. Finally, a rotation 
by an angle <|>2 about takes the B-based coordinates to the general 

position shown by ^ in Figure 9. 

Now, the matrix 0 transforms vectors based In the Inertlally fixed 
coordinates to the B-based coordinates. So 6 Is given by: 
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cos iJ'j cos 0 


cos >J )2 0 -sin 


0 

0 

1 

e = 

-sin i/'j cos ili^ 0 


0 10 


0 cos sin 


1 

0 

0 

1 


sin ® ’('2 


0 -sin cos 


’ cosi)/ 2C08 i)J 2 (cosilJj^slniJ/j + simJfj^slniJi^cosiJ^j) 


-C08l|i28tn<l)2 (costtij^coaiJ)^ " 8lmlij8lni))2Bliu()2) 


Blntli- 


-8lmpj^cos\J)2 


(B±n'l>j^8±arli^ 

-COBtpj^Blnip^^OBtllj) 

(slnilijCosxlJj 

+ cosl(/j^sliu(i28lnil<2) 

COSlfl^COSlJlj 

( 1 . 94 ) 


The rotation coordinates of the body A are defined relative to the body 
B, In Figure 10 . Flrat, a rotation by an angle 6^ about b^ brings a^ 

In the position shown by F^. Then, a rotation by an angle 62 about ?2 
brings the axes a^ In the position shown by G^. Finally, a rotation by 
an angle 6^ about brings the A-based coordinates In the general posi- 

tion shown by a^ In Figure 10 . 

Let 6.„ be the matrix that transforms vectors In the A-based coordl- 
AB 

nates to B-based coordinates. Then 0 ^ Is also given by Eq. ( 1 . 94 ), but 
In which the angles <li2 and are replaced by 8^, 62 and 8^. 

Now, taking 8j^ and 82 small but keeping 8^ arbitrary, 0 ^ Is given 
by 




0088^ 

slnSj 

-02Cos8 

-sln82 

cose^ 

8j^cos8 

®2 

-«L 

1 
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site to those between the bodies B and C. He now refer to Figure 8 for the 
following analysis. 


Let and be the contact force and torque vectors, respec- 

tively, applied by the body C on the body A, and are expressed In the A- 

based coordinates. Let ^ and ^ be the force and torque vectors ap- 

“*BC 

piled by the body C on the body B and are expressed In the B-based coordinates. 

Let and be the shift of the center of mass of the body B In B-based 

coordinates and that of the body A In A-based coordinates. The vectors 

and defined In Figure 8, are expressed In A-based and B-based 

coordinates, respectively. In this section, the terms sbc. 

will mean the vectors F^_, C expressed in the B-based coordinates. 

~Av# “Titi A 

Let 6j^, 02 and 9 ^ be the relative angular positions of the body A 
with respect to the body B. It Is assumed that 6^ and 62 are small, and 


6^ Is large. It Is also assumed that the slopes of the body C at the end 
of A with respect to the end B along the axes b^ and ^2 are given by 62 
and respectively. 


be expressed as 
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Then 


and 


where 




T 


" lAi.C 


~ " ’^BA^b^AC^ 


(1.98) 


(1.99) 


If s* T * r 

b^A * 

From Eqs. (1.96) and (1.98), we have 

he = - + ^3®J - 1^7 W W + ^9®^ * 

Now let ll^ll^ total masses of the bodies A and B. Then 

because of the force-free environment of the satellite, we have 

+ IlMlfeCg = 0 . 


Therefore 


where 


jAi ^ ~ y-^ 


IWIa 

y = 

ll«l^ 

Therefore, from Eqs. (1.100) and (1.101), 


(1.101) 


( 1 . 102 ) 


Similarly 


(1.103) 


he “ ^®AB^^1 " ■*■ ^®AbW ” ^8^®AB^^ 


+ te;^(T7 - rre)SABl^A + ^ ®I^ ^9®> 


1®IJW - Y^2>®Ab’^ + t®I^<^7 - ^"8^®AbJ^ + (®^ ^3® + ®aS ^9i> 

(1.104) 


.-1 


o-X 
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In the above equation, 0^ transforms vectors In the A-based coordinates 
to vectors In the B-based coordinates. The neglected term Is due to small 
angular accelerations. From Eqs. (1.96), (1.99) and (1.101), 


T 

-=BC 




^ Y ^'^10 ^ ^'^11 * ’^ ab '^ b ^^^ 


- ^AB^l’i - t^l2+ 'aB^ 


(1.105) 


and 


4c = ^ - y^ii>®ab14v 


■*■ ^®AB^’^10 ■ ^’^11^®AB^^ ®AB^’^6^ ■*■ '^12®^ 

— ^®AB^"^4 “ ^5^®AB^^ ■*■ ^®AB^^10 " ^’'^11^®AB^4v ®AB^’'^6® ''' ’'^12®^* 

(1.106) 

Kow let the body C be taken to be a uniform shaft having E^, and 1^, 
as Its modulus of elasticity and the area moment of Inertia of the cross- 
section, respectively. Let 2^ be the length of the shaft. It Is now 
assumed that the displacements of the shaft at the ends A and B are equal 
to and C^, respectively. It Is also assumed that the rotations of 

the shaft at the ends A and B, about the axes ]b^ and b^, are the same 
as those of the bodies A and B, respectively. Then 
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where A^, 


The symbol 


[Tl] 






^C®C 


( 1 . 107 ) 


r-tj] 


iTl] 


( 1 . 108 ) 


[Tsl 


6E I 
C C 


~P~ 

X 


( 1 . 109 ) 


rj = nt^. Tg = TIT2 


and 


T, = TIT3 


( 1 . 110 ) 


Is the cross-sectional area of the shaft, and 

* 


= i" ' 


( 1 . 111 ) 


stands for the modulus of viscosity of the shaft material. 


[T4] 




!Vc 

0 


( 1 . 112 ) 


[Ts] 


[T4] 


( 1 . 113 ) 
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0 


0 




[T6] 


^Vc 


(1.114) 


[Tiq] = nlT^l, [T^j] = nCTj], [Tj2^ = n[Tg) + 0(^3] . d-us) 

Substituting the expression for 6^ and the Eqs. (1.107) through (1.115) 
in Eqs. (1.103) through (1.106), the values of the contact forces and 
torques are obtained. 


10. The Assembled Equations 

In this Fart I, the equations of motion of the different structural ele- 
ments have been obtained. The derivations are shown for the body B. A simi- 
lar set of equations can be obtained for the body A. 

The configuration shown in Figure 1 has been described by 110 general- 
ized time-dependent position coordinates. Accordingly, a total of 110 
equations are generated for the body B. This set of equations in matrix 
notation has the form 

+ [A 2 (<D,ai,t)]{qg} + [As((>),w,t)]{q^} = {Si,(oi),ai,t)} . (1.116) 

Similarly, the equations for the body A take the form 

[AjKq*} + [Aj(ci),a),t)]{q*} + [Aj(w,w,t)]{q*} = {Aj(u,a),t)} . (1.117) 

In these equations, [Sj] and [aJ] are (110 x 110) square matrices of con- 
stants. These final forms of the equations will be used in the next part 
of this work. 
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11. Error Bounds of the Solutions 


The complete analysis of this report Is based primarily on the assumed 
solutions of the partial differential equations for the beams, plates and 
shells. These solutions obviously are not exact. The Galerkln's method 
only minimizes the error corresponding to the assumed form of the solution. 
A scheme Is now devised to estimate this error based on the assumed solu- 
tion. 


In each of the cases of beams, plates and shells, the solution u Is 
sought for an equation of the form 

Au = pu + f (1.118) 

where A Is a linear differential operator, p and f are functions Inde- 
pendent of u. Let V be the error In the assumed solution, which Is 
given by (u + v) . Now, let (u + v) actually satisfy an equation 

A(u + v) = p(u + v) + f - Aj^(u + v) (1.119) 

where A^ also Is a linear operator. Then from Eqs. (1.118) and (1.119), 
we obtain 


or 


Av = pv - Aj^ (u + v) , 


(p - A)v = Aj^(u + v), 
or 

V = (I - P ^A) ^p ^Aj(u + v). 

Now, It Is assumed that p ^A Is a contraction operator. Then 

V ^ (I + Ap ^)p ^Aj^(u + v) 
or 

||v|| = 11(1 + p"^A)p"\(u + v)l| . (1.120) 

In the Eq. (1.120), Aj^(u + v) la the quantity obtained by substituting 
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the assumed solution In the governing differential equation. Also A and 
Vt ^ are known operators. Then the maximum values of the error v Is ob- 
tained by taking suitable norms In the Eq. (1.120), over the complete time 
and spatial domain of the operators A and y. A very easily calculated 
norm is In the space L^, and this norm will be used. Hence, the follow- 
ing operations are required: 

a) Carry out the complete analysis, and obtain the coefficients of 
the Galerkln functions, 

b) evaluate the norms In Eq. (1.120), 

c) If ||v|| Is too high. Increase the number of Galerkin's functions 
and repeat the procedure, 

d) If ||v|{ Is small, then the assumed solution Is satisfactory. 
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Nomenclature 


(Note: Unless otherwise mentioned, all vectors are In B-based coordi- 

nates. For symbols not defined In this list, see "Nomenclature" In Fart I.) 




^s^lB^ 




[L3],[LJ 

M 

B 

e, Bj 

^BM 

^BC 


Uo 

* 


J 

a 

[Ib] 

’■Bl’ ^B2’ ^B3 
^EBC,1 W=1.2.3) 


Continuity matrices, defined In Eq. (2.3) 

Continuity matrix, defined In Eq. (2.12) 

Continuity matrices, defined In Eq. (2.20) 

Continuity matrix, defined In Eq. (2.48) 

Continuity matrices, defined In Eqs. (2.73) and (2.73a) 
Residual magnetic moment vector 
Element of defined In Eq. (2.74) 

The Earth's magnetic field vector 
Elements of defined In Eq. (2.75) 

Magnetic torque vector on the body B 

Eddy current torque vector on the body B 

Permeability of the materials In the body B 
Velocity of light In vacuum 

Position vector of an element of the body B from the 
center of mass 

Volume eddy-current density vector of the body B 
Static electrical conductivity of the body B 
Moment of Inertia matrix of the body B 

Diagonal elements of [I ] 

Elements of TgBC’ Eqs. (2.81-2.83) 
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p 


T 

-4ibs 


P 

e 

eo 

n 


jr 

aoi t)o, CO 

^EBG 

d 

<ii. di. da 
U 


R 




* 


Tg^. i = 1,2,3 


^B 

[Pg,]. i = 1-4 


4 


•> Average density of the elements of the body B 
o Solar radiation torque on the body B, In B-based 

coordinates 

-7 2 

= 1 X 10 Ibs/f t for a surface normal to the Sun 

° Reflection coefficient of the body B 
= Unit outward normal vector to a surface element of 
the body B 

° Unit vector directed from the Sun 
° Elements of at t = 0 

= Gravity gradient torque vector on the body B 
= Unit vector towards the center of the Earth 
= Elements of defined In Eq. (2.92) 

° The Earth's gravitational constant, 1.4082 x 10‘ 
ft’/sec* 

° The distance between the center of the Earth and 
the body B 

= Angular momentum vector of the body B 
= Position vector of a mass element of the body B 
= Angular momentum vector of rigid rotors or reaction 
wheels Inside the body B 
= Control torque vector on the body B 
* 

= Elements of defined In Eq. (2.102) 

— D 

<=■ Total torque vector on the body B 
° Matrices and vector, defined by Eqs. (2.104) and 
(2.105) for the body B 

= Angular momentum vector of the body A In A-based 
coordinates 
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* 

*-A 




T 

-£A 


[P^ll. 1=1-4 
[Aj^]. 1 = 1-4 
[A^] ,1 = 1-4 
iMg^]. 1 = 1-5 
[M^], 1 = 1-5 


= Angular momentum vector of rigid rotors and reaction 
wheels Inside the body A, In A-based coordinates 

= Control torque vector on the body A In A-based 
coordinates 

* 

= Elements of T. 

—A 

= Environmental torque vector on the body A In A-based 
coordinates 

= Total torque vector on the body A In A-based 
coordinates 

= Matrices and vector for the body A corresponding to 
[P , ], defined by Eq. (2.110) 

Bi 

= Matrices and vector for the body B, defined by 
Eq. (2.112) 

= Matrices and vector for the body A defined by 
Eq. (2.113) 

= Matrices and vectors for the body B defined by 
Eqs. (2.115a) and (2.116b) 

= Matrices and vectors for the body A corresponding 
to defined by Eqs. (2.123) and (2.124). 
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1. Introduction 

We have obtained only the equations of motion of the Individual struc- 
tural elements In Part I. Now we shall use the continuity of displacements, 
rotations and moments at the junctions to reduce the dimensions of the prob- 
lem. By eliminating the generalized position coordinates, equations are ob- 
tained Involving only the angular velocities. 

Equations (1.1) through (1.120), referred to in this Fart II, were 
given in Part I. 

2. Continuity Conditions 

In Che analysis thus far, only a few of the continuity conditions be- 
tween the various elements of the structure have been considered. The re- 
maining conditions will be considered now. 


(a) End Mass - Beam Junctions 

The condition that the sum of Che beam-end mass inertia forces and the 
beam-end shear forces must be zero has already been Introduced In the deri- 
vation of the equation of motion of these masses. The remaining conditions to 
be Imposed are that the deflections of the tip-masses must be equal to the 
beam-tip deflections, l.e., 

^ •fb3j}s=, (2.1) 

for the (l,j) pairs given by (17,1), (18,2), (19,3) and (20,4). For ex- 
ample, let the case 1 = 17, J ° 1 be considered now. 

Here [vi®]^ Is the Identity matrix. So for 1 = 17, j <» 1, Eq, (2.1) 
becomes 


'n,i 

“ ’1.1 


'17,2 

“ l0l,2(®‘’**-^^ 

s'*!, 2 

'17.3 


2‘*1,3* 3**1,3 
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Thus, for each beam there Is a relation 

where [,L ] are constant matrices. More details of this equation for 

D Ic 

this and other end masses are given in Appendix 11. 

(b) Rigid Body - Plate Junctions 

The force and moment continuity conditions between the rigid bodies 
and the plates have already been considered in Eqs. (1.81) and (1.87). The 
displacement and rotation continuity equations are now obtained. 

The displacements (y^, 1 = 1,2, 3, 4) of the centers of mass of the 
rigid bodies are taken to be the average of the plate displacements at the 
four points at which the rigid bodies are fixed to the plates. For the 
plate function which is singular at a Junction point, the function is eval- 
uated at a small distance e away from that point. As an example, let 
the conditions for the rigid body of mass m^ be considered here. 

The four Junction points of rigid body No. 1 and the plate No. 2 have 
their position coordinates in the plane of the plate given by aj, as, 
as, and aio. Then, referring to Fig. 7 and Eqs. (1.35) and (1.40), we 
have 

yi = ^IX 2 (a 7 ) + X2(“8) + X 2 (as) + X2(oii))] (2.4) 

with 

ttj ■= tj exp(10j) . (2.5) 

Therefore 
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yi 




+ (rio-l)rio] + j^g[(r9-l)r| + (rio-l)rio] 

+ a^ coa 09 + Ji(Xirio) cosBio] + g[(r9~l)r9 cos 09 

+ (rio-l)rio cos0jo] + a^ gt(r9-l)ri cos 09 + (ni)-l)rio cos0io] 

+ ^2 j^[J2(X2r9)co8 209 + J 2 (X 2 rio)cos 20io] + ^2 cos 209 

+ (rio-l)rio cos 20io] + a^ j[(r9-l)r| cos 209 + (rio-l)rfj cos 20jo]} 


+ 


1 

4 


10 

I 


k=7 


X "2.J '^2.j(V • 


( 2 . 6 ) 


Similarly, 

y2 = X 2 j^(t) '*’^^2 

+ ^2 j^^[(rii“l)r 1 J + (rj2-l)ri2] + ^2 + (ri2"l)ri2] 

+ ^2 2lJi(Xjrii)cos0ii + Ji(Xiri2)cos0i2] 

+ &2 g[(rii"l)rnco 80 n + (rx2-l)ri2 cos0i2] 

+ 82 g[(rii-l)rii CO80M + (ri2-l)rf2 cos0i2] 

+ 82 j^[J2(X2rji)2 CO80I1 + J2(X2ri2)2 cos0i2] 

+ &2 2l(rii-l)xii COB 2011 + (ri2-l)ri2 cos 20 12) 

+ 82 2((rii-l)rfi cos 20n + (ri2-l)r|2 cos 2612]} 


1 14 18 

+ j I 5^ X 2 j *2 • 

k=ll j =7 *■' 


(2.7) 
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The rotational continuity conditions are now obtained for the rigid 
body No. 1. Again, referring to Figure 7 and Eqs. (1.3S), (1.40) and (1.77), 
we have 

®1 1 “ ■ X2(a?) - X2(a8)l (2.8) 

* 0 

®1 2 “ + X2(“0 - X2(oi 7) - X2(aio)] . (2.9) 

Similarly, for the second rigid body, we have 


2,1 2i, 


■[Xs(“ll) + X2(“12) - X2(«I3) - X2(0H4>] 


6 


and 


®2 2 “ 2t“IX2(au) + X2(“>t) - X2(“i2> - X2(“i3>] 


( 2 . 10 ) 


( 2 . 11 ) 


Combining Eqs. (2.6) through (2.11), the displacement and rotation continuity 
conditions for the two rigid bodies are obtained as 
* 

(rSl p<>2 (2.12) 



where [ L„] is a matrix of constants, 
r B 


(c) Beam-shell Junctions 

The In-plane rotation of the shell and therefore the Inplane moment 
have been assumed to be zero. The force continuity in the radial direction 
has already been considered while deriving the equations for the beam-end 
masses. The continuity conditions for the shear forces in the transverse 
direction of the beams are neglected, as the shell displacements in those 
directions are taken to be zero. The remaining continuity conditions for 
each beam-shell junction are three displacements, two rotations and two 
moments . 
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(i) Displacement relations — Referring to Eqs. (1.3) and (1.59), the 
displacement relations at the junction point and 6 = 0^ are 

given by 

(2.13) 


3'*1,2 “ ° 
3**1,3 “ ° 


(2.14) 


^2^1 . . ,...P3Pl, 


ll.l ” (*1+P')(zi-P")[bj,(t)e ^ ^ + b3(t)e 

ir(P'+P.) 2ir(p'+p ) 

+ b^(t) sin — j COS0J + b 3 (t) sin j cos6^ 

TfCp'+P^^) 2ir(p'+p.) 

+ bg(t) sin j sln6^ + b^(t) sin ^ sln6^ 

ir(p'+p.) 2ir(p'+p ) 

+ bg(t) sin j cos 26 + bj(t) sin j cos 29^ 


ir(p'+p^) 


2ti(p'+pJ 


+ bjg(t) sin sin 20^ + b 3 ^(t) sin j — — sin 26^ . 


(2.15) 

(11) Rotation relations — The conditions that the slopes of the beams 
at the joints are equal to the slopes of the shell there, are 


a. r 1 Zlail 

2'*1,2 ’’l‘l^l,2^ ■ a- 30 

B J*l-Pi 


e-6. 


1 ®^r"l 

2’l,3 *’l^l‘‘l,3^ ■ ~3^ ^ 


from idilch 


e-0. 
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2*^1, 2 


2‘*1,3 


1 ®^Br‘ 
■ a ~36~ 

O 


r 1 1 ^^B.r 


Zl=pi 

0=6^ 


"l“Pl . 
6=6. 


(2.16) 


(111) Moment relations — The bending moments on the i*'*' beam at the 
junction are given by 


and 


“2 “ “ ^b®l^^b^l,3^Pl ^1^1,2^ 


M 3 = - 3)p3(j^qj^^3) 


(2.17) 


The bending momenta on the shell at the Junction In the tangential and 
axial directions are given by 


= 


b^B • s^B ^.'^B,r , 

ft n • 


3 " 5 . 




2^_2 " s>'B ^2 


12(1 - sV^B 


3zt 


1 -jVPi 

0 = 0 . 


and 


sS • sS ^%,r 


+ .U. 


2-t„2 a, 2 8 B ggz 


12(1 - s ^ B ^^ ^==1 


_lVPl 
6 = 0 . 


(2.18) 


Now, setting 


[“2] 

L J 

r«el 

1 1 


H f 

1 

1 

1 M 


(2.19) 


the moment relations are obtained. 
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Combining the Eqs. (2.13) through (2.19), the complete shell-beam 
junction continuity conditions for the 1*^ beam are given by the relation 

“ fs^lBlfs<»B> ^2.20) 

where constant matrices. 

(d) Plate-shell Junctions 

In this problem the force continuity conditions at the plate-shell 
Junction are neglected because the plate and shell displacements are assumed 
to be zero there. Only the displacement, slope and moment conditions will 
be set up. The corresponding equations will be obtained for each harmonic. 

(1) Displacement relations — Let the Eqs. (1.34) and (1.59) be used to 
obtain the relations for the junction between the plate No. 1 and the 
shell B. 

As the In-plane plate displacement Is zero, so 


Again as the shell axial displacement is zero, so 

= 0 . ( 2 . 22 ) 

As all the other terms of the Eqs. (1.34) and (1.59) are zero on the bound- 
ary, so the shell and plate displacements are matched at the junction. 
Similarly for the boundary condition at the other end of the shell, one must 


X2,i(0 “ 0 • (2.23) 

(11) Rotational relations — At the boundary. It Is assumed that the 
joint Is rigid, and the plate and the shell rotate by equal amounts. There- 


fore the relations are 


1 _ 1 ^^B.r 

‘ “b **^Jr-l “ “b ®*iJ z,-p" 
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To evaluate the Eqs. (2.24) and (2.25), the following result is to be noted; 


M 

3 r , 
J r«l 


2(1 - a^a^) 


where 4 'j ^ (t,0iOi,) is the plate displacement function exemplified 

1,J J 

by Eq. (1.38). 

Equating the terms Independent of 6 in Eqs. (2.24) and (2.25), we ob- 
tain the relations 


-Xoa^^^3J,(Xo) + a3_3^ + a^^^j + + 2 (1 - 

= -)i.i[b2(t)eP*‘^" + b3(t)eP’P"] 


-A(ia2,i3Jiao) + * ®2,15 * ^ " 'P^2,j 

= 21 ( 82 ( 1 ) 6 '^=“^’ + b 3 ( t ) e “ P * P ’] • 


Similarly, equating the coefficients of cos6 In Eqs. (2.24) and ' 
(2.25), we obtain 

a^ ®i 8 ®1 9 “ (2.29 


33 yXiJc(\i) + ^2 8 ^ ^2 9 ” ” ' 


Equating the coefficients of sln6, we get 


nil Jo (A i) + a, ,, + 


1,11 “ 1,12 2 i 


^ lb « - 2 b 7 ] 


32 , 1 oAiJo(Ai) + a^ + 33^32 “ ‘ 
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The coefficients of coa 26 are equated to give 
*1 ^ 2 “l 3 “ (2.33) 

®2 ®2 2 ^ ®2 3 ^ ~ ^ * (2.34) 

Similarly, the coefficients of sin 26 generate the relations 
*1 '*’ ^1 5 ^ 6 ~ 2bji] (2.35) 

®2 ^ ®2 5 ®2 6 “ ” '*’ • (2.36) 


Sow from the Eqs. (2.21) and (2.27) through (2.36), the shell displace- 
ments are given by 


bi(t) = 0 


(2.37) 


b*(t) = i[e<P®P"-P®^’> - e(P®P"'P®P’>3'^ jo-P®'^’{-X.Jt(Xo)a^^^3 

^ 11 

®1,14 * ®1,15 ■*■ ^1,7 ^ Jg ~ 

„ 18 ) 
+ eP®P (-loJi(l,)a2^33 + + 2 (1 - rj)X2,3 3| 


(2.38) 


b3(t) = i [e^P^P"-P®P'^ - e<P®f’"-P®P’^r^/e-P®P'[-XoJ.(Ao)a^^ 


.13 


11 


■*■ “l,14 *1,15 \,7 ■*■ ^ 

n" 

+ e**® [-AoJi(ls)a2^33 + ^2,U * ®2,15 * ^ 


b-.(t) = - »z,7> + =1.8 + \,9 ~ “2,8 “ “2.9^ 


(2.39) 


(2.40) 
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86 


bs(t) = - + ®2,7^ “ 1.8 “1,9 ®2,8 “ 2 . 9 ’ 

be(t) = - aj.lO^ ■*■ “l,ll ''' “l,12 " “2,11 ” “2,12^ 

(2.42) 


b9(t) 


‘2,10> + “1, 

,11 * 

“1,12 ■'■ “2,11 

“2,12^ 




(2.43) 

> “1,2 + 

“1.3 

■ “2,2 ' “2,3^ 

(2.44) 


■^“1. 

3 ■*■ “2,2 ■ “2,6^ 

(2.45) 

? * “1,5 

“1.6 

" “2,5 ' “2,6^ 

(2.46) 

,4> + “1,5 “1,6 

“2,5 “2,6^ 

(2.47) 


Equations (2.37) through (2.47) la rewritten as 
where 


„ *, * * T 

tp** ^ • p'’2> 


and “ ‘^°"®bant matrix. 


(ill) Moment relations — The moment relations are obtained by summing 
to zeto the radial bending moment on. the plates and the axial bending mo- 
ment on the shell at the shell-plate junctions. Thus the relations become 

1 1 , 3\~ 

I| *p“l ^■*'p‘'l^r 3r ^Jr=l 

E • r 3*5., 

s B SB ^B.r , ^ ^B.r 

12(l-^y|)a| 3z! “ ® 38* zi=p" 
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L 


—r + 

3r^ P ^ 


.1 ^ ^ 

' r* 38^ _ r=l 

s^B • r ^ ^^^B,; 

!(1- u|)a| L 3*1 ° ® 30^ Jzi=-P’ 




To evaluate Eqs. (2.50) and (2.51), the following results are to be 

used. If ijj (r,6,a) are the plate displacement functions, one of which 
1*3 

Is shown In Eq. (1.38), Chen 

3\ i"| , (1 - r.^)^ 

= 2(1 - rj) + 4 J (2.52) 

3r* Jr=l J 1 - 2r^co8(e-6j) + r* 


TT 1 - r: 


-TT 1 - 2rjCO8(6-0^) + 

V (1 - r|)cose 

-IT 1 - 2rjCos(6-6j) + r* 

V (1 - r*)sln6 

-n 1 - 2rjCos(8-9j) + r* 

I T (1 - r*)co8 29 

-TT 1 - 2rjCos(6-0j) + r* 

, IT (1 - r?)sln 20 


d0 “ 1 for 0 < r . < 1 


d0 = TjCOsOj 


d0 = r^slnOj 


d8 “ r? cos 20. 


-TT 1 - 2rjCos(6-6j) + 


d9 = rf sin 20. 


2(1 - r*) 


2(1 -rj). 
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The values of Xj^t X2 ^ ^ given by Eqs. (1.34), (1.40) and (1.59) 

are now substituted In Eqs. (2.50) and (2.51). Then the tetms independent 
of 6 and the coefficients of different harmonics are equated separately. 
Equating the terms Independent of 6 In Eq. (2.50), we obtain 

p“lW,13^^"pV^'’'^*^’''’^ “l,14<^+pV ^1,15^%V Xi^y(3+pP^) 

11 


+ I [2(l-r?)(3+2 u )x, .]} 

j*8 ^ P A ±,j 


= - — 2^[2(l4p2fi)eP^P"b2 + 2a4p3fi)eP'P"b9] 

12(1-3^1) 

Similarly, from Eq. (2.51), we obtain 
pD2{a2,i3(l-py2)>^oJx(Xo) + ®2,14^%V ^2,15^®V2^ 

18 


(2.59) 


^ [2(l-rj)(3 + 2py2)X2,jD 


s B s B 


[2(1 - p2li)e‘P^P’b2 + 2(1 - p3fi)e"P®P b3] • (2.60) 


To obtain the equations for the harmonic components In Eq. (2.50) and 
(2.51), these two equations are multiplied successively by cos6, sln6, 
cos 20 and sin 20 and Integrated with respect to 0 between the limits 
(-Tr,ir). This operation on Eq. (2.50) generates the following four equations: 


p“l{^>Jo(^i)(pV‘l - + P^'1^^1,8 + p’^1^^1,9 


+ 8 [(1 - rj^)rjcos0j Xi,j 


]} 


(2.61) 
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pD^{A»Jo(Xi)(pV^-l)a^^^O + (2+py^)aj j + (i+pVa^ 

11 

+8 I [(1 - r?)r 8in0 X, .]> - 0 

Vl‘Dai.1 ■" <Vl^*1.2 ^VP“l,3 

11 

+8 I [(1 - r?)r cos 20 v, <]} “ 0 

j„8 J J J ■‘■•J 

pDj{X2Ji(X2)(pV^-l)a^^^ + (4+pUi>a^ 5 + (6+pyi)a^ 5 
11 

+ 8 I [(1 - r?)r Bin 20 x, J) = 0 . 

j„8 J J J ■‘■•J 

In a similar way, the following four equations are generated from 
Eq. (2.51): 

k|Jo(Xi)(py2 - D^2,7 * pV®l,8 + + P^2^“l,9 

18 


(2.62) 


(2.63) 


(2.64) 




(2.65) 


XiJo(Xi)(py2 - 1)82 10 + + p»'2>“2.11 + p'^2^^2.12 

4 8 [(1 - >;)r; .IrfjXjjl - 0 

X2Ji(X2)(pU2 - 1)82 1 + + p‘'2>°2,2 + + p’"2^“2.3 


(2.66) 


18 

+ 8 I [(1 - rj)rj cos 20 jX2^j] 


J“7 


(2.67) 
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A2Ji(X2)(py2 - Da2,4 P^'2^^2.5 p^2^®2.6 


18 

+ 8 t(l - r^)rj^ sin 20 jX 2 ^j] = 0 . 


( 2 . 68 ) 


This completes the derivation of the continuity conditions for the individual 
elements. 


(e) Remarks on Reduction of Dimensions 

We shall determine the effect of the continuity conditions on the total 
number of the degrees of freedom in our complex structure. 

From Eq. (2.48), it is seen that the shell coordinates can be expressed 


as a linear combination of the coordinates of the two plates, i.e.. 


(2.48) 


(2.69) 


From Eqs. (2.20) and (2.48), we have 

^b‘«*i> = ls^iBJ^s4> = • 

This shows that the beam coordinates can be expressed in terms of the plate 
coordinates. The quatlon (2.12) states that the rigid body coordinates can 
also be expressed in terms of the plate coordinates, which are 


* * T r 

fr’l* r‘»2^ = to; q } 


(2.70) 


From Eqs. (2.3) and (2.69) the coordinates of the beam-end masses can be ex- 
pressed in terms of the plate coordinates by the relation 




(2.71) 


All these relations show that only the spring-mass-damper coordinates 
cannot be expressed in terms of the plate coordinates. 

When Eqs. (2.38) and (2.39) are substituted into Eqs. (2.59) and (2.60), 
we obtain two relations between the plate coordinates themselves. So 
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Eqs. (2.59) through (2.68) and Eqs. (2.22) and (2.23) supply 12 dependence 
relations between the 28 -t- 21 = 49 coordinates for the two plates. This 
means that only 49 - 12 = 37 plate coordinates are Independent, and the rest 
are dependent. These dependence relations are defined by an eqtiatlon of the 
form 

= tpS’^p<»> (2.72) 

where {^q} Is a (37 x 1) column and ® ^ constant matrix. 

Since there are four spring-mass-damper coordinates, the motion of the body B 

Is defined by the (41 x 1) vector {q_l. So the originally mentioned 

B 

(110 X 1 ) vector {q_} Is given by 
B 




[Lb] {93} 


(2.73) 


where [L^] Is a (110 x 41 ) constant matrix and Is composed of matrices 

W' ‘a^lbJ* tb^lBl W- 


The equation for the body A Is similarly given by 




'■•.IV 


(2.73a) 


3. External Torques 

It has been assumed that the external forces on the satellite are negli- 
gible. The center of mass of the satellite remains static. But the 
existence of external torques cannot be neglected. The basic method of 
modelling the external torques follows that shown by Dobrotln, et al.[l] 
and Tidwell [2]. These torque models are derived for the body B. The 
torques for the body A are derived similarly. 

(a) Residual Magnetic Torque 

The flow of electric current In the electric circuits Inside the space- 
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craft reacts with the Earth's magnetic field to create this torque on the 
spacecraft. 

Let ^ be the residual magnetic moment vector of the body B In the B- 
based coordinates. Let B be the Earth's megnetlc field vector In the B- 
based coordinates. The elements of ^ and ^ are assumed to be of the form 

{Mg} = [M, M, M]''' (2.74) 

{b} “ [e cos (Dg gt, -e sin tOg ^t, B^]^ • (2.75) 

Then the external magnetic torque on the body B Is given by 


T 

-EBM 


i^B 


X B 



(B^ + e sin 0 )g ^t) 

(e cos ojg^gt - B 3 
-e(sln oig jt + cos Wg ^t) 


(2.76) 


(b) Eddy Current Torque 

Eddy currents are Induced In a satellite spinning in the Earth's mag- 
netic field. These eddy currents in turn react with the surrounding mag- 
netic field to create a torque on the satellite. The eddy current torques 
on the body B are given by a volume integral over the body B, as follows: 


T 

-EBC 


/ r x(j X B)dv 

Use 


(2.77) 


where 

Uo = the permeability of structural material of the satellite, 
c* = the speed of light in vacuum, 

r^ = the position vector of a volume element from the center of mass 
J = the volume eddy current density. 

All vectors in Eq. (2.77) are In the B-based coordinates. The vector ^ Is, 
in turn, given by 
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(u X B ) X r + V<ti 


(2.78) 


— 2lloO '^B 

where a >= the static electrical conductivity of the satellite material, 

^ ° the operator for spatial gradient, 
and i)i = a scalar potential for the body such that V*<() = 0 and = 0 
on the bounding surfaces of the satellite. 

In this problem the field of la the thin plates ahd shells. Therefore, 

Ip Is nearly a constant. This reduces Eq. (2.77) to the following form; 


T 

^BC 


2p?a. 


/ r x{[(ux^)xr]x B}dv 




2yoOc 


2yo^c 


T 

•%BC 


/(B • ‘ • 


(2.79) 


2bo'Jc 


It Is now assinned that the moment of Inertia matrix of the body B has the 
form 

“^Bl ° ° 

0 I 


fv = 


B2 

0 0 1 . 


0 

B3 -■ 


(2.80) 


Hence Eq. (2.79) reduces to 
1 


-EBC,1 - 2y§oc% ^“b,i"^“ ^,3^ + 2 “b, 2®^" 2a>3_3t)I. 

+ (Bimg^3 - B 30 .g^ 3 e cosu.3^ 31)133] 


B2 


(2.81) 


189 





=BC,2 2vhc*p 

+ (b 1 Oig_2 + 8 ^“b,3^ «^“"b,3'^)Ib.31 


(2.82) 


+ (e\,3sln^ + B,e<a3^2«ina.3^3t)Ig2] 


where 


T 

"H;BC 


fT T T I 

^ EBC.l’ EBC,2* EBC,3' 


(2.83) 


(2.84) 

and P la the average density of the structural materials. Using the condi- 
tions that <»)g and <o^ 2 are small compared to 0 )^ 3 and e Is small 
compatsd to B 3 , the Eqa. (2.81), (2.82) and (2.83) are simplified to 

f “b. 1>^B,3] 


T 


<;-(«’«^B.3«l”V3t B, m3^2)V3j' 


(2.85) 


(c) Solar Radiation Torque 

The solar radiation torque arises from the asymmetric pressure distribution 

developed on the surfaces of the satellite due to the electromagnetic radiation 

from the Sun. The formulae used here are taken from Beletskll [3]. 

The radiation torque vector, T^„_, on the body B, Is given by the fol- 

^tSoS 

lowing Integrals over the surface of the body B exposed to the Sun. 

^gg ■= Pg{(l - eo)[l + eo [2 / n x r^(n- 2 )*ds]} ( 2 . 86 ) 

In Eq. (2.86), all vectors are In the B-based coordinates, and 
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-7 2 

p s a constant = 1 ^ 10 Ibs/ft for a surface normal to the Sun, 
e 

to = the reflection coefficient of the body, 
a = the unit outward normal to the surface element, 

T the unit vector directed from the Sun, 
r^ = the position vector of the surface element from the center 
of the body B. 

In this analysis, it is assumed that the angle of Incidence of the solar 
radiation Is large on the plates, which are also highly reflective. This 
assumption allows us to consider that the plate surfaces produce negligible 
torque, so the Integrals In Eq. (2.86) are taken over the surface of the 
shell only. 

Let 

T 

To “ (ao, bo, col 

be the direction cosines of the sun-vector at the time t = 0. Then 



COS U)_ 

B, J 

sin < 0 ^ 3 ! 

0 



ao 

T *» 

-sin <c^_ 3 t 

cos (^ 3 ! 

0 



bo 


0 

0 

1 

_ 


_co_ 


(ao cos ojg ^t + bjj sin lo^ ^t) 
jr = (-so sin oig ^t + bo cos ^t) • 

CO 


(2.87) 


Let the unit outward normal on the shell element be given by 
n = [co 80 , sln 0 , 0 ]^ . 


Then 


T = co80(aocos 


B.3 


T 0| 




B.3‘ 
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Now, using Eq. (2.88) and taking 

T 

^ = [agCoaS, agBlne, BgZj] 

Eq. (2.86) Is reduced to the following form: 


^BS = 


tao[(p") -(p')*]-liCoH)p ,t}{-aoslnWg ^t + boCOsUp ,t} 


B.3 


B,3‘ 


{liCo«)g^3t - ao((p") -(p’)*]HaocosUg gt + boslntOg ^t} 


- -f Pg-ea-a^[(p")^'(p')^l(2a? -i- b§) 


B.3" 

B.3" 


(2.89) 


(d) Gravity-gradient Torques 

The gravity-gradient torque vector on the body B, In the B-based 


coordinates. Is given by 


Iebg = f, td- l,d] 


where 


- “ td^.d^.d^] 


the unit vector towards the center of the Earth, 
the Earth's gravitational constant 
1*4082 X 10® ft^/sec^. 


and 


R the distance of Che center of the Earth from the body B. 
Then, from Eqs. (2.80) and (2,90), we have 


(2.90) 


-EBG 




^’•BS ' ^B2^‘^2*^3 
^^Bl ' Sl^®1^3 
^^B2 ' ^Bl^‘^l*^2 


(2.91) 
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Now, let the Inertlally fixed reference coordinate axes be set up such 
that Xj points In the direction perpendicular to the orbital plane and along 
the angular momentum vector, X^ Is along the line joining the centers of 
the masses, and X^ Is along the orbital velocity vector. Then In the X- 
frame, located on the body B, ^ Is given by 
^ = [0, 1, O]"^ . 

Therefore In the B-based coordinates, ^ Is given by 
d = 0.d^ 

where 6 Is the transformation matrix given by Eq. (1.94). 

Carrying out this transformation, we obtain 


dj^ = cosiJij^slnij/j + slnijij^slnilJjCosiiij 

= cosif/j^costlij - slnij/j^slntji^slnijij (2.92) 

dj ° - slnij)j^co8i(/2 • 


Substituting the Eqs. (2,92) In Eq. (2.91), and using small angle approxi- 
mations for and <|< 2 > we obtain 

f<^B3’^B2>t-^^2®^“'''3 - (’''l + ¥A - 1 


T 

-EBG 




(lB2-lBl^ti^^'^'^l)®^" 2 i(- 3 + i(-^t|-2COs2i))3] 


(2.93) 


(e) Total Environmental Torque 

In this analysis only the above-mentioned environmental torques are con- 
sidered. Then the total external torque, T„„, on the body B Is obtained 
as 


T »T +T +T +T 

■4:b -ebm -ebc -4:bs -ebg 


(2.94) 
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These torques do not Involve the {q.} as the effects of the flexibility on 

O 

them are Ignored. 


4. Moment Equations for the Composite Bodies 

The moment equations for the body B are derived now. The corresponding 
equations for the body A can be obtained In an identical manner. 


(a) Angular Momentum 

Let ^ be the angular momentum about the center of mass of the body B 


In B-based coordinates. Then 

N 




(2.95) 


where r^ Is the position vector of a differential element of mass dm, ex- 
pressed In the B-based coordinates, from the actual center of mass of the 
body B. Since 


r 

-c 




(2.96) 


we obtain 




/(Cg + r) [(^+£)+<^’‘ (c^+r)]dm 


where r^ is the B-based position vector of the element from the nominal 
center of mass, or 


+ • (2.97) 

Because c„ and c^ are assumed to be small, and since 
/r dm = - ||M||gC^ 

Eq. (2.97) Is linearized by keeping only the fourth and the last terms. The 
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last term can be vritten aa h^} • In this expression is 

the angular momenta of rigid rotors or reaction vheels inside the body B. 

Thus Eq. (2.97) reduces to 

Therefore 

or 

\ “ ^b“b ■*■ ^b“b ^ 

+ “gljlOg + + 0^ X ( /r X £ dm) . (2.99) 

Neglecting 1 , and using the symbols from earlier parts of this analysis, 

O 

we get 

■ Vb + “bVb ^b “b*^ * I 

Ji^bPi I [(£^)nil<*«i + b^A I 

+ I pPfct //p\ Xk dS + // p\ XfcdS] 

+ sPb'// ^bV® + // (2.100) 


where 

J » 1,2 and 13-20, k = 1,2. 

Introducing the expressions for i]^, ^ and ^ obtained before, 

Eq. (2.100) is transformed into 


Ifli-B + “B^B!i!B + hg + + iPg^llg + 


(2.101) 


Details of Eq. (2.101) are shown in Appendix 11. 
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(b) Equations of Motion of the Body B 

* 

Let the controlling torque vector T on the body B be expressed as 

fi 


* , * * ft ^ 

Tg = Tj2* ’’^83^ 


( 2 . 102 ) 


In B-based coordinates. Then the total torque applied on the body B In 

“D 

B-based coordinates Is given by 


h “ Ib+Ieb+Ibc+ 


(2.102) 


In which shown In Figure 8. In this equation, 

Is obtained from Eq. (1.103), T_„ from Eq. (1.105) and from Eq. 

”BL> “TLJJ 

(2.94). Therefore 

h = 4 + Ieb + <"4 + ^AB^1> - ("5 + ^Ab"2> 

+ "l - "2 )’Vb% <\0 + ^AB^7> - (\l ^^AB^8> 


^'^12 '^AB^‘^9^- 


(2.103) 


The equation of rotational motion for the body B Is then given by 

•"d 

^B “ ifMb • 

Hence from Eqs. (2.101) and (2.103), we get 
^*’b1^% ^“b*’^B1^ ■■ ^AB^7^ ” ^^11 ■*■ ’^AB^B^ 

+ {[(^5 + ^^^ 2 ) “ y(T4 + ’^ab”'^!^ “ ®BC^Y 

= .4 ^EB - ^ ^ ^AB>^3li - I"l2 + + ^AB>^9^® 


- - *1b ' “b^!b - %h% (2.104) 
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Equation (2.104) is rewritten as 

^ ^ f^3l% = ^*’b4^ ' <2.105) 

In this equation Is a (3 x 1) vector. The symbols 

[Pgg] stand for (3 x 41) matrices. The elements of tPgj] <^B3^ 

constants, but those of [P„-] involve o)„ . 

dZ — fi 


(c) Equations of Motion of the Body A 

N" 

Let ^ be the Inertial time derivative of the angular momentum vector 
of the body A In the A-based coordinates. Then the equation corresponding to 
Eq. (2.101) Is given by 


4 “ Va + Va“a “a*’a '^ai^’a “a^^ai’^b • 


(2.106) 


Let ^ be the controlling torque vector on the body A In the A-based 


coordinates, so that 
* 

^A “ 


* * * , 

<'^A1’ ’^A2* ’^Al^ • 


(2.107) 


Also, let Tg^ and ^ be the external and the total torque, respectively, 

on the body A In A-based coordinates. Then, referring to Figure 8, we have 
* 


T = T+T +T +(R xp). 
-A ^ '^C ■^C'' 


Substituting Eqs. (1.104) and (1.106) In the above equation, we get 
A -1 

h = t®AB<\ - ^5>®AB 

^AC®AB<’^1 " ’''^’^2^®AB^®A^A% <®AB<^10 ‘ ^^11^®AB 

®AC®AB<’^7 " "'^8^®AB^®A^A^ <®AB<’^6^ '^12^^ 

+ ^AC®I^<"3i-^^4>J <2-108) 
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(2.109) 


Then the moment equations for the body A are given by 
—A 

Hence, from Eqs. (2.106) and (2.10H), the Eq. (2.109) reduces to a form 
given by 

^ “ fV ' (2.110) 

More details of Eqs. (2.10S) and (2.110) are in Appendix 11. 


5. Equations for the Angular Velocities 


(a) A Review of Various Equations 

The behavior of the complex dynamic system shown in Figure 1 has been 
formulated into six sets of equations. These equations are: 



h 

(1.116) 

- • ••* - , •* - - * 

^1 ^A + ^2 “Ja *3 \ 


(1.117) 

^8 “ ^**8 


(2.73) 

4 “ ^a'»a 


(2.73a) 

^81^*8 ^82‘*B *'b3'‘b 

*"84 

(2.105) 

^Al^A *’a2^A *’a3'*A 

” ^A4 

(2.110) 


These equations are linear in the generalized position coordinates 
and 

analysis so far has been along well-known techniques, but now a completely 
new approach will be taken. 


but nonlinear In the angular velocities ^ and The 


(b) Condensation of the Appendage Equations 

The general appendage equations (1.116) and (1.117) Involve very large 
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matrices. To try to solve these equations would require a great computa- 
tional effort. The usual practice [9], [10], [11] Is to truncate the co- 
ordinate vector by a considerable amount. Instead of gross truncation, 
the Eqs. (2.73) and (2.73a) are used now to reduce the dimensions of Eqs. (1.116) 
and (1.117) In the following analysis. 

Equation (2.73) Is first substituted Into Eq. (1.116) to obtain 

° ^4 * (2-lU) 

Now If Is an (m x 1) vector and If Is a (n x 1) vector, with 

m > n, then the Eq. (2.111) Is a system of m equations In only n Inde- 
pendent variables. The number of equations Is reduced by premultiplying 
Eq. (2.111) by (L_)^. Thus we obtain 

A 

(2.111a) 

This equation Is rewritten as 

^l^B *B2‘*B ^3 ’b “ *B4 * (2.112) 

Similarly, using Eq. (2.73a), Eq. (1.117) Is reduced to 


(2.113) 


(c) Derivation of Angular Velocity Equations 

He now eliminate ^ from Eqs. (2.105) and (2.112). It Is to be noted 

that, as P . are rectangular matrices, their inverses do not exist. Hence 
Sx 

P of Eq. (2.105) are augmented by adding (n-3) rows of zeroes and then 

DX 

added to the square matrices In Eq. (2.112) to form the equation 

(Aji + Pgi)qg + (Ag2 + ^b2^^B ^^3 ^B3^‘*B “ ^*B4 ^B4^ ' (2.114) 
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The matrix + Pgj^) is square with definite inverse. 

Equation (2.114) is now pre-multiplied by ^ and sub- 

tracted from Eq. (2.112) to get 

^^B2 ■ ^l^^l *’b 1^ ^^2 “ ^l^^l ■*■ *’b1^ ^^3 ■*■ ^B3^^‘*B 

^■^4 ” ^Bl^^l ^Bl^ ^^4 '*’ ^B4^^ ' (2.115) 

The above equation is rewritten as 

«Bl’B ^ “b2<1b = «B3 • (2.115a) 

Differentiating Eq. (2.115a), we obtain 

«B1^B + K + «B2^^B ^ •Sl-^B = S3 * 

In this equation the value of q„ is substituted from Eq. (2.112) to get 

D 

”bi(^1^4 ~ ■^B1^2^B " SiSs’b^ ■*■ %1 ■*■ *S2^^B ^2‘*B “ ^3 

or 

(^1 “b2 ~ “b1*Bi\2^‘*B * (^2 ■ *Sl^B1^3^'*B ° ^^3 " “b1^1^4^ * 

Here q„ is replaced by its value from Eq. (2.115a), and we obtain 
B 

(^2 - «Bl^k3 - (^1 “b 2 - «BlSk2>«;K2J'lB 

“ (^3 - “b ^4 - V «B2 - »B • (2.116) 


Equation (2.116) is rewritten as 


^4'^B “b5 

<’b = ^4X5 • 


(2.116a) 


(2.116b) 
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Substituting the value of q fiom Eq. (2.112) In Eq. (2.119), we get 

15 

f'Sl ^2 “ “b1*B1^2 ■ “b1*B1^3 “ ^^1 ^2 " “b1^iS2^^1^2^’b 

■'■ ^^2 " ^lSl^3 “ ^^1 ■*■ ^2 ■ “b1*B1^2^SiS3^'’b 


f^3 “b1^1^4 ^^1 ■'■ “b 2 " “b1^1^2^Si*B4^ ' 


( 2 . 120 ) 


Now, substituting Eqs. (2.115a) and (2.116b) Into Eq. (2.120) and simplify- 
ing, we obtain 

^^Bi“b2^^3 ■*■ ^^^2 " ^3“b^1 ■'■ SA^1*B1^3^“b2^'^3 

+ t^3"^2n3 = "S4 • (2.121) 

In deriving Eq. (2.121), It has been assumed that [^2^ small compared 
to and [A^j], and so only the terms linear In A^^Ag 2 snd ^^ 3^52 

were retained. On further simplification, Eq. (2.121) reduces to 
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“ ^B3*B3^^2*^2'^3 ■*" “b3 
^^Bl'^Bl ■ *’b 3^3^S4 • 


The equatien (2.122) Is the haslc differential equation In the angular ve- 
locities only. It la a general equation where the effects of any number 
of flexible elements can be taken Into account. Beginning with the equa- 
tions (2.113) and (2.110), a similar set of equations for the body A can also 


be obtained In the forms 


and 


= m.Hi 


A4 AS 


(2.123) 


- "a3\>aA 2’«A3 ^ t<^A2^l2 - ^3^I^>^2“I^“a3 ^ «a3 


^^Al^Al * ^A3^A3^^A4 


(2.124) 


In the next part of this work we will begin with these highly nonlinear 
equations and look for special cases and their solutions In transient. 
Intermediate and long-time ranges. 
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M2 “ 0; F2A’^ > P^A^^ 


b) Case 2: 0; < Pj^A"! 
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i) 

Case 

3 - 1 ; 

M2 ^ 0; 

P3A-/ > 
-1 

P A'^- 
-1 

Mj^ = 0 

-1 -1 

il) 

Case 

3 - 2 ; 

M2 / 0; 

-1 

PlA, ; 
-1 

W " "1^1 
-1 -1 

ill) 

Case 

3 - 3 ; 

M2 ^ 0; 

P3A3 > 

PlAi ; 

^2^2 " ^"^1 


d) Solutions for the special cases 
Stability Criteria of M^ 


Conclusions 
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Nomenclature 


i - 1,2,3 
i = 1-4 
P^, i = 1-4 

a 


1 

2 

(ij 


*ij’ ®lj 


d 

n = 0,1,2,... 

h 

2. 

u 

V 

f 


P.Q 



= Same as or of tlie prevloua report 

= Same as or of the previous report 

= Same as P^^ or P^^^ of the previous report 

° Independent set of the generalized position 

coordinates (see Part II) 

= 4 

= Identity matrix 

“ [<i^, or [ 0 ^, 

= The largest element in (Aj^a”^) 

= Matrices, defined by Eqs. (3.6) and (3.7), 
respectively 

= Vector defined in Eq. (3.34) 

= Components of w, defined in Eq. (3.11) 

= Component of (O, defined in Eq. (3.12) 

= ^ 

= 

= [d,0]'^ 

= Pu + qd 

= Matrices, defined in Eqs. (3.30) and (3.31), 
respectively 

= Operator, defined by Eqs. (3.32) and (3.33) 

= Solutions of u 
= Cartesian vector space x i 
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6 


k 


M 


T 

* 

01 




V 

* 

i (s.z). 


T 


= Inf [u(t) - ^(0) I for 0 ^ t ^ T > 0. 

= A constant, defined In Eqs. (3.36) and (3.39) 
= A constant, defined In Eq. (3. 37) 

= Outer limit for t 

= Component of u, defined by Eq. (3.43) 

= Component of defined by Eq. (3.44) 

= Gradient operator with respect to ^ 

“ Gradient operator with respect to u 

= Vectors defined by Eqs. (3.52) and (3.53), 
respectively 

° Enlarged time scale, defined in Eqs. (3.56) 
and (3.67) 
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1. Introduction 


In this part of our work we begin to explore the consequences of Eqs. 
(2.122) and (2.124) derived earlier. These are very general and highly non- 
linear equations. For convenience, Eq. (2.122) together with Eqs. (2.112) 
and (2.115a) are repeated here without the subscript B as follows: 

- P3A"^)A^m"^]M3 + - P3A"^)A2M2^]M3 + M3 = (P^a"^ - PjA"^)^ 

(2.122), (3.1) 

Aj^a + Aji + A3£ “ ^ (2.112), (3.2) 

«xi + “ 2 ^ “ ^3 (2. 115a), (3. 3) 

First, we shall obtain perturbation series solutions for the generalized 
forces M^ and valid for different time zones. Then we shall con- 

sider some special cases when the structural parameters combine In a way to 
simplify these equations to a great extent. 

In what follows, the subscripts A and B are dropped from the symbols 

“a1* ^1’ *A1’ ^1’ *’a1’ *’b 1’ % analysis Is similar for 

either of the two bodies A and B. 

2. Asymptotic Solutions of M 3 
a) The Equation 

We now look for the asymptotic solutions of the general eqtiatlons (2.122) 
or (2.124) represented by Eq. (3.1). For the present, we assume that the 
generalized stiffness matrix A 3 Is large. Noting that Involves 

A 3 ^, Eq. (3.1) can be represented as 

+ eB^j(<£.t)yj + y^^ = edj^(u,t) (3.4) 

where 
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i = 

M3(u,t) 



e = 

the largest element In (Aj^a”^) , 


^13 

e ^[(P^Aj^^ 

- P3A3^)Aj,m"^] 

(3.6) 

^3 = 

e“^(P2A’^ 


(3.7) 

0l = 

e-\PiA-l - 


(3.8) 

0} = 

[i4. 

or [ojg, i^] . 

(3.9) 


This shovs that remain quantities of comparable magnitude even 

though e Is small. The form of the Eq. (3.4) makes It clear that as the 
structure becomes Increasingly rigid, e tends to zero and the Eq. (3.4) 
approaches the equation 


y^ = 0 . (3.10) 

on further simplification, Eq. (3.10) becomes the rigid body equation 
for the two bodies A and B. 


b) The perturbation series 

We begin the process of solving Eq. (3.4) by assuming (s In the 
following form: 


2 3 

( 0 ^ = + e + e + ... (3.11) 


Defining 


we get 


hf = + .. 


“l “Oi **i • 


(3.12) 


(3.13) 


Since 


1 = = ^(“q + h, t) 
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We now expand 2 , 1''^ * Taylor’s scries sround 
Defining 




'1.3 3w, 


and 


we obtain 


i 


’'l.jb “ 3Wj8iSj. 




( 3 , 14 ) 

Therefore 




( 3 . 15 ) 


Bence 


jf^(w,t) « yi<«o.t) + e(yjj(«o.t)»ij + 2yi,j%»'=)% 

+ yi,j%’*>Sj K.jk%’">“ij^k ■*- 2yi,3k^j“ik ^ yi,jk%’*^%3^k^ + 

( 3 , 16 ) 

Expanding A^(w,t) in a Taylor’s series about Mq, we obtain 
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or 


Similarly, 

Bij(i!i.t) = B^j(^.t) + 

We also have 

dj^Cu.t) = d^(^,t) + ed^ j(^,t)iDy 

■*■ ^ ^‘*1,1 ■*■ I ‘*l,jk%’*^^“lj“lk^ ■*■••• 


(3.17) 


(3.18) 


(3.19) 


The series given by Eqs. (3.14) through (3.19) will now be used to 
generate the perturbed equations from which u will be solved. With our 
choice of e. It Is seen that each expansion series approaches the corre- 
sponding nominal value as the satellite becomes Increasingly rigid. This 
characteristic Is in agreement with the physical nature of our problem. 


c) The perturbed equations 

Substituting Eqs. (3.14) through (3.19) in Eq. (3.4), and separating 
0 12 

the coefficients of e ,e ,e we obtain the following equations: The co- 
efficient of gives 

y^(^,t) = 0 . (3.20) 

The coefficients of e give 

or 

y^j (“o’*^^“lj “ • (3.21) 
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From the coefficients of c , we get 

+ [yi.j(<i!o*‘)“2j + 1 yi.jk(iSo*‘>“ij“iki 

“ ■*■ 2 ‘*l,jk^“0‘*^^“lj“lk 


■'■ 2 yi,Jk^^'*^^“lj“lk “ ‘‘"2 ^1, Jk^^'*^^“lj“lk 


2«i,Jk - yi,jk>%’‘^“lj“lk 

Now, differentiating Eq. (3.21) with reapect to u, we have 

•"i.) “ yi,kj“ik + ^i.j • (3.2 

Differentiating Eq. (3.23) with respect to u, we obtain 

^l,jk ■ ^l.jk “ 5'i^jki“u ■*■ • (3.2 

Substituting Eqs. (3.23) and (3.24) In Eq. (3.21), and neglecting 
dj(t^,t) and dj (j^,t), we get 

^l,kj(^’*^^“lk“2j ■'■ 2 yi,jkl(^’*^^“lj“lk“lf ''' 2 ^1, jk(^'*^^“lj“lk ° ° 

Then neglecting the cubic term In the above equation, 

“2j = "Kj • (3-2 


Then from Eqs. (3.11), (3.21) and (3.25), 
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(3.26) 


“l “ “Oi 

where ^ Is to be obtained from Eq. (3.20). 


3. Convergence Criteria of the Asymptotic Solutions of 


a) Selection of the time domain 

In view of the simplifications made In solving Eq. (3.4), It may be 
argued that the obtained solution of that equation Is quite different from 
the actual solution. To check If such Is the case, we proceed In the fol- 
lowing way: 

Let 

yj^ = • (3.27) 

Then Eq. (3.4) 1s given by 

£A(u,t)x = ed - eB(u,t)x - . (3.28) 

Let the vectors ^ and v be defined by 

u = [x.il^ 

and 

V = [d,0]'^ . 

Then Eqs. (3.27) and (3.28) are represented by 

u = Pij + = f(^(t),t) (3.29) 

where P and Q are given by 


and 




0 


0 


(3.30) 


(3.31) 
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(3.32) 


From Eq. (3.29), the formal solution for ^ Is 

t 

u(t) = u(0) + f(u(t),t)dt 

or 

u(t) = F[u(t)l (3.33) 

vfaere the operator F Is defined by the Eqs. (3.32) and (3.33). 

Now our problem can be stated as follows: 

Suppose u^(t) Is the exact solution of Eq. (3.32) and 0^(1) Is an 
approximate solution, such that both solutions have the same Initial value, 

l.e. 

u(0) «> Uj^(O) = U2(0). 

Then If [F(iij^) - ^(112)] Is small, will [u^(t) - U2(t)] be also small? 

If the answer Is In the affirmative, then we also want to know the condi- 
tions under which It Is so. It can be proved by applying the Contraction 
Mapping theorem [1], that [U2(t) - U2(t)] Is small If F Is a contraction 
operator. We now obtain the conditions for F to be a contraction operator 
In our problem. 


Let IR be the 2n-dlmenslonal Cartesian space containing the n- 
dlmenslonal vectors x and Let ^ belong to a set U such that 

u 6 u = {u : u e K, |u(t) - u(0) | £ 6} (3.34) 

for 0 £ t £ T > 0 

and 

||u|| = max |u(t) I (3.35) 

te[0,T] 

where 5 Is a positive constant. 

Mow It Is assumed that there exists a constant k such that 

|f(u^j^,t) - f(u^2ft)| - “2I (3.36) 
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where 


— 1’— 2 ^ ^ ^ ® 10,T]. 

It Is also assumed that there exists a positive constant M such 

|f(u,t)| < M, u 6 U. t € [O.T] . 

Now, from Eqs. (3.30), (3.31) and (3.32), we have 

iKij^.t) - f(u2,t)| = |PUj^ - Pu^ + Qd(Uj^,t) - Qd(u2,t) 


- “2^ " ^ 2^1 

< |P + Q*Vdl-|u^ - U 2 I 

A"^(Vd-B) - i a"^ 

“ ^ • l“i - “ 2 1 ' 

I 0 ^ 

Then comparing Eqs. (3.36) and (3.38), we get the constant k as 


k «■ the greatest eigenvalue of 


A"^(Vd-B) 

I 



From Eqs. (3.32), (3.3A) and (3.37) we have 
t 

|u(t) - u(0)| £ |f(£,t)|dt £ MT . 

From Eqs. (3.32), (3.33) and (3.36), we have 

t 

||p(^, ) - F(u )|| = max (/ [f(u„t) - f(u,,t)]dt| 

^ ^ te(0,T] ” ^ 

t 

< max / k|u. (t) -u,(t)|dt 

t€[0,T] “ ^ ^ 

< kT llu^ - U2II . 

Then from Eqs. (3.40) and (3.41), we see that the operation 


that 

(3.37) 


(3.38) 

(3.39) 

(3.40) 


(3.41) 

F will be 
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a contraction if we choose T such that MT < 5, kT < 1. Therefore 

T<mln{|,i} • (3.42) 

Thus, £q. (3.42) gives the time interval in which our approximate solutions 
will remain close to the exact solution. This also shows that the system 
is unstable if k is very large. 

b) Truncation of the solution series 

An Important piece of Information we require in this analysis is the 
number of terms of the series given in £q. (3.11) that are necessary for 
the solution to be arbitrarily close to the exact one when t f [0,T]. We 
now show how to estimate any term of that series from the value of the pre- 
ceding sequence. 

* 

Let the partial sums of the series in Eq. (3.11) be denoted by 
such chat 

u = y e“u (t) . (3.43) 

n n 

n»0 

it 

Let be defined by 

Then Eq. (3.4) can be expressed as 

or 

Cl = + ^Ci ■ C ' 

where is the gradient operator with respect to or 

Cl = eliC-t) + Vd(^,t)(Ci - C - 

where ^ is the gradient operator with respect to m or hr . 
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Therefore 


[I - = ed(^,t) 

- e[^(^,t) + A(^,t) + B(^,t) -^] ^ . 

— n ^^2 flt 

with n = 0, Eq. (3. AS) becomes 

[I - e^(^,t)]i* = e^Co^.t) 

or 

= e[I + e^(^,t)]d(^,t) 

provided eVd(eiQ,t) ^ -j . Similarly 

[I - e7d(M*,t)]y* = ed(u*,t) - e[7d(M*,t) 

, * , d* * d * 

+ A(Wj^,t) ^ Zi . 

Substituting Eq. (3.A6) in Eq. (3.47) and dropping the 
we get 

[I - e^(w*,t)]i* = ed(u*,t) - e^(^(i£*,t) 

+ A(w*,t) ^ + B(w*,t) ^)d(i!^,t) 
or 

i2 = + eW(u*,t)]d(u*,t) 

since 

4(<4.0 = 0 

80 

d(u*,t) = • 

Therefore 

^2 = + eVd(ai*.t)][7d((u*.t)]y* . 


(3.45) 

(3.46) 

(3.47) 

term involving 


(3.48) 
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Carrying on this procedure we get 

y*) ]^ = e[I + eVd(ai*.t)HVd(tii* ]^.t)](y* - y * j^) for n^l . (3.49) 

This shows that after choosing T, we have to choose n such that the 

* * 

difference of ~ obtained from Eq. (3.49) becomes small for 

t £ [O.T]. 

c) Validity of asymptotic solutions - Jump conditions 

We now consider Eqs. (3.27) and (3.28) to obtain another Important 

characteristic of that system. These two equations are represented by 

. * 

ex “ f. (x,.z) 

i “ (3.50) 

where e Is a small parameter. The property that e appears with the 
highest derivative of the state vector makes It a singularly perturbed 
system. In the previous section we have assumed that the system should tend 
to a definite limiting mode as e approaches zero. We will now examine 
whether or not an asymptotic solution of the system exists. 

To understand the problem, we consider the phase velocity vector ^ 
given by 

V = = [-^ l*(x,x), B,(x,x)l^ 

or 

V = f*(x,i),0]'^ + [0. £(x.i)]^ • (3.51) 

The second vector on the right-hand side of Eq. (3.51) does not depend on 

* 

e. The first vector becomes infinite as £ approaches zero if ^ r 0. 

This means that when the phase point P(x,x) Is not on the surface 

* . 

1. ^ the component ^ of the phase velocity away from the surface 

la great. At the same time, the component ^ is limited. Such motion Is 
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maintained until P(jc,^) comes very close to the surface f (x,^) = 0. 
After that the variables change with finite velocities. However, for some 
the point ^7 again move away from this surface and the system 

may lose equilibrium. This phenomenon makes the singularly perturbed 
systems difficult for analysis. 

The equilibrium state of the system is given by the degenerate system 
f*(x,^) = 0 (3.52) 

£(x.i) = Z * (3.53) 

If the system is at equilibrium, then k can be solved as 

X = x(i) (3.54) 

from Eq. (3.52), and then Eq. (3.54) can be substituted into Eq. (3.53) to 
yield an explicit solution for 

But if the system is not at equilibrium, then it will not be possible 
to obtain Eq. (3.54) from Eq. (3.52). This means that the jump condition 
is given by the relation ^ 

31 

det [-5^ ] » 0 . (3.55) 

Using this relation in our problem, we see from Eq. (3.27) 

L (i»Z) “ 6A - Bx] - A • 

Hence * 

3f 

dim [-5 — ] = 0 for all t . 

e -»■ 0 - 

But though Eq. (3.55) is satisfied, our system is in equilibrium, as 

11m X = 11m f^ (iiZ) “ ~ - — A 

e-»0 e-*-0^ e-*-0 ^ 

= A ^(d - B^) < 0= 

since 
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11m ^ = 0 . 

e ->■ 0 

So It Is seen that the Jump condition does not exist In this problem, and 
the asymptotic solution Is applicable for all time. 


A. The Inner Boundary Layer Solution of 

We now consider the governing Eq. (3.4) of our system and Its basic 
solution given by Eqs. (3.20), (3.21) and (3.25). It Is evident that the 
asymptotic series has reduced the differential equations to algebraic re- 
lations. It has happened because the time scale chosen for these solutions 
is unity. This means that these solutions are valid for large values of 
time. So, the time scale Is now enlarged near t = 0. The boundary layer 
solutions thus obtained exhibit marked changes in their nature. 

Let T be given by 

t = ex for 0 ^ t ^ e • (3.56) 

Then, expanding A^j(u, t) around t = 0, and w = we get 

+ eA^j(<£,0)T + -Y Aj^j(u,0)x^ + ... (3.57) 


A^jKO) = A^j(^,0) + eAy^^(^,0)u^^ 




= A^j(c^,0) + £A^j^^(^,0)a)^^+ ... 


Similarly, B^j(<o, t) and d^(u, t) are expanded in Taylor's series around 
t = 0 and u = ^ • 

Then, substituting these expansions In Eq. (3. A), and using the 
relations 
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dt 


(3.60) 


e dT 


and 




df 






dT 


(3.61) 


the coefficients of e ^ and e° generate the following equations: 

— Z(Un,T) = 0 (3.62) 

dT^ 

[A(7y) 1(^.0) + B(o)q,0) Ip x(“0*0> <^.63) 

dT 


Hence, for T < 1, we have 




i.(<^,T) = aT + b (3.64) 

and 

co^(T) ■= -[(V i)‘^ A'^](o^.O)[B(UQ,0)a + b]T (3.65) 

where a,b are constants and 

Uj^(O) = 0 . 

Thus at the inner boundary layer , u is given by 


u = ^ + (3.66) 

where is the solution of Eq. (3.64). The constants a and b will 

be obtained by matching this solution with the asymptotic solution at 
T = 1 . 


5. The Outer Boundary Layer Solution 

A necessary condition to be satisfied by all boundary layer solutions 
is that these solutions must coincide with the large time solutions as t 
becomes infinite. It is seen that the basic large time solution for ^(<£,t) 
is given by 

i(^,t) = 0 . (3.20) 
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The corresponding Inner boundary layer solution Is given by 

« at + b . (3.64) 

As these two solutions are Incompatible as T we look for an outer 

boundary layer solution with a different time scale. 


Let T be given by 

t ° 

Tci for 0 ^ t ^ I'T » 

(3.67) 

With this scale, we have 


1 

(3.68) 

dt 



1 ^ 

(3.69) 

dt* 

e dx* 



The functions A(u,t), B(m,t), ^(<o,t) and ^(w,t) we expanded as before 

around u ° ^ and t ° 0. These expansions and the relations (3.68) and 
(3.69) are then substituted Into Eq. (3.4). Then, collecting the terms con- 
taining e° and /e, and equating that to zero, we get 

Aij(<!^,0)yj(<^,T) + 71 B^j(^,0)yj(^,T) + yj(uQ,T) = 0. (3.70) 

Equating the coefficients of e to zero, we get 

[A(^ i)] (^,0)ffij^ + 2[A(^ ^)]l£j + V[a£ + 2,](^,0)u^ = ^((^,0) (3.71) 

Thus It is seen that with this time scale, the quantities i(ji)g,T) and 
^j^(^,T) satisfy the equations for damped oscillations. Setting x = " 

In these solutions leads to the long term solutions. Equations (3.70) and 
(3.71) are solved by using u(0) Instead of ^ In the coefficients of 
yy yj. “i and . 

6. Special Cases 

We now consider some special cases of Eq. (3.1). 
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a) Case 1: M2 = 0 

If for a given set of parameters, the nominal value of M2 Is either 
zero or very small compared to Mj^, then from Eq. (3.3), we have 

M^i = M 3 (3.72) 

Or 

i = ^^^3 . (3.73) 


Differentiating Eq. (3.72) and neglecting M^, we get 
q = . 

Integrating Eq. (3.73) by parts, and neglecting we have 



(3.74) 


(3.75) 


Substituting q, q, and q from Eqs. (3.74), (3.73) and (3.75) into 
Eq. (3.2), we get 

AiM ‘^3 + A2M'^3 + A3M'^ M3dt ■= A^ . (3.76) 


We now consider the three possible variations of M^ when M2 is zero. 


1) Case 1-1 ; M2 = Mj^ = 0. In this case the only solution for 
(3.1) from Eq. (3.3) Is 

M3 = 0 . (3.77) 

Again, 


«1 = M2 = M3 = 0 

Implies 


0 • 


This means that the angular velocities are given by the rigid body motion 
of Che satellite, and are not Influenced by the flexibilities of the 
elements . 
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In this case 


II) Case 1-2 : ” 0 and ^2^2^ '' 

(I + '*' ^ ^ Identity matrix. Therefore 

" A"^[I - (I + P2A^^)~^(I + 

= Aj^tl + (I + P^a“^)“^(I + P 2A"^) 

+ (1 + PjA^^)"^(I + P2A2^)^ + .•■] 
or 

^ 3a"^[I - P^A~^ + • (3.78) 

Then Eq. (3.76) reduces to 

A^A'^d - P^A^^ + P2A“^]M3 + (I - P^a"^ + PjA'^IMj 

+ AjA'^U - P^a'^ + P 2A'^] M3dt = y . (3.79) 

Now, usually [A2] is small. Therefore, Eq. (3.76) becomes 

(M3A'^A3M'^)N3 - (M^a”^A^) (3.80) 

vhere 

t 

N, = / M,dt . (3.81) 

3 *^0 3 

Then using Eq. <3.78), Eq. (3.80) becomes 

N3 + 3[(P2A"^ - P2A“^)A2A'^A3A“^]H3 = 3(P3A"^ - P2A”^)A2A'^A^ . (3.82) 

III) Case 1-3 ; Mj = 0 and P2 A“^ > P^a”^ . In this case, 

(I + P2A“^)'^(I + P2A2^) > I . 

Then 

= -A'^d + P2A2^)‘^(I + P3A"^)[I - (I + P2A2 ^)"^(I + P^a"^)]"^ 
or 
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(3.83) 


S -A2^[21 - 3P2A"^ + 3P^a'^] , 
Then Eq. (3.76) reduces to 


Aj^a“^(2I - 3P2A"^ + 3Pj^A”^)M 3 + (2l - 3P2A“^ + 3 Pj^a“^)M3 
1 1 1 

+ A3A2 (21 - 3P2A2 + 3P3Aj^-^) M 3dt = -A^ . (3.84) 

If [A2] Is small, then Eq. (3.80) becomes 

N3 + 2[(P2A“^ - P3A]^^)A2A"^A3A2^]N3 = (P^A^^ - P2A2^)A2A"^A^ . 

(3.85) 


b) Case 2: M2 ^ 0 and P3A“^ < P^A^^ 

If M2 0, then the general equation (3.1) is to be considered. A 
simplified form of M2 Is derived now. We have 
(I + P3 A"^) < (I + P^a"^). 

Therefore 

= a'^U - (I + P 3 A'^)"^(I + PjA"^)]"^ 

= a"^[I + (1 + P^a'^)”^(I + P 3A"^) + (I + P3A‘^)“^(I + ^34"^)^ + ...] 

or 

M2 ^ S Sa'^EI - P^a'^ + (3.86) 

1) Case 2-1 : M2 0, P3A”^ < P^A'^, M^ = 0. From Eq. (3.3), we 

have 

M2q = M3 (3.87) 

or 

q = m ”^3 • (3.87a) 

Differentiating Eq. (3.15) and neglecting M2, we get 

M24 = *3 (3.88) 
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and 

M2q ■= ■ (3.89) 

Substituting q, q, q from Eqs. (3.87), (3.88) and (3.89) In Eq. (3.2), we 
get 

= A^ . (3.90) 

Then from Eqs. (3.86) and (3.90), for small A2, we get 

M3 + 3 M 2 A"^(I - Pj^A^^ + P3A2^)M3 = M2A“^A^ 

or 

M3 + - P3A3^)A3Ai^^“3 “ ~ 

II) Case 2-2 ; M2 * 0; P^a"^ < P^a"^; P 2 A"^ < P^a"^. 

The Eq. (3.86) Is substituted In Eq. (3.1) to obtain 

3[(P3A"^ - P3A3^)A3A"^1M3 + 3[(P2A‘^ - P2A3^)A2A'^]M3 + M3 = (P^a"^ - P 3A'^)A^ . 

(3.92) 

This Is the required equation for In this case. 

III) Case 2-3 ; M2 0; P3 A^^ < P^a"^; P 2 A“^ > P^a"^. Here the 

form of the reduced equation does not change and Eq. (3.92) remains valid. 

c) Case 3: M2 ^ 0; P3A'^ > P^A^^. 

In this case, 

(I + P3A'^) > (I + P^A^^) • 

Therefore 

= a"^[I - (I + P 3A^^)“^(I + P3A3^)]‘^ 

= -A3^(I + P3A2^)"^(I + P^A"^)!! + (I + P3A3 ^)"^(I + P^a”^) 

+ (I + P3A3)”^(I + P^A^^)^ + ...] 
or 
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(3.93) 


= -Aj^rZl + 3 PjA^^ - 3P3A"^] , 

I ) Case 3-1 : ^ 0; P 3 A 3 ^ > Pj^A^^; = 0. For small k^. 

Substituting Eq. (3.93) In Eq. (3.90), the equation for Is obtained as 

M 3 - M 2 A^^( 2 I + 3P3A^^ - 3P3A“^)M3 = 

or 

M 3 + [2(P3A"^ - P 3 A^^)A 3 A^^]M 3 - (PjA”^ - P 3 A‘^)A 3 A"\. (3.94) 

II) Case 3-2 ; Mj |t 0; P 3 A“^ > P^A^^: < P^a"^. The Eq. 

(3.93) Is substituted In Eq. (3.1) to obtain 

[2(p3A“^ - PiA"^)AiA^^]M3 + (2(P3A“^ - P2A"^)A2A3^]M3 + M3 

= (P^A'l - P3A"^)A^ . (3.95) 

This Is the required equation for M 3 In this case. 

III) Case 3-3 ; M 2 ^ 0; P 3 A"^ > P^a'^; P 2 A"^ > P^^^. In this 
case, too, the equation for ^ remains the same as Eq. (3.95). 

d) Solutions for the special cases 

The long-time solution and the boundary layer solutions have been ob- 
tained for Eq. (3. 4 ) which Is the general equation for the system. By 
redefining the variables A^^ , and the equations for all the 

special cases derived earlier can be brought In the form of Eq. (3.29). 
After that the analysis for the special cases Is Identical to that of the 
general case. It Is Interesting to observe that the basic solutions for 
oiQ(t) are not affected by those simplifications except In the cases 1-2 
and 1-3. 
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7 . Stability Criteria of 


In our present notations, the Eqs. ( 2 . 122 ) and ( 2 . 124 ) are given by 


[(P^A^^ - P3A3^)A^M"^]M3 + 1(P2A2^ - P3A3^)A2M"^]M3 + M3 
= (P^a"^ - P3A‘^)A^ . 


( 3 . 1 ) 


It is to be noted here that the motion is stable when the angular velocities 

^ and ^ are stable. Now M3 la a complex differential form Involving 

01 . and So for the motion to be stable, a necessary but not sufficient 

"A 

condition Is that M3 should be bounded. 

Examining Eq. ( 3 . 1 ), It Is evident that as P^, P3, Aj^, A3 and A^ 

are bounded, M3 will be bounded If and only If the coefficients of M3 and 

M3 are positive semldeflnlte. It Is also seen that A^ Is positive definite 
as It Is the generalized mass matrix. So this makes the coefficient of M3 
posltlve-semldeflnlte. Then the remaining required condition is that 
((PjAj^ - P3A3^)A2M2^J must be posltlve-semldeflnlte. 

How, M2 Is positive If ^ ^2 negative if 

Pj^Aj^^ < P3A3^. Considering our previous analysis, we see that If P2 and 
A2 are basically negative while A^, A3, P^ and P3 are positive, the 
coefficient of M3 will be positive-definite If the following Inequalities 
hold: 

PiA"^ < P3A3^ < P2A'^ ( 3 . 96 ) 

PjA"^ > P3A“^ > P2A“^ . ( 3 . 97 ) 

If 


M 


1 


0 , 


then from Eq. ( 3 . 18 ) we see that stability is ensured If A2 Is posltlve- 
semldeflnlte. 
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Again, If 


= 0 , 

then from Eq. (3.4), the stability is also ensured if is positive- 

semldeflnlte. 

If all the quantities Pj^, Pj. **3* Aj and A3 are positive, then 
the stability criterion for is given by the following inequalities: 

Pj^a"^ > PjA"^ > P 2A"^ (3.98) 

and 

P^a"^ < P 3A"^ > . (3.99) 

8. Conclusions 

In this report, we have analyzed several characteristics of the system 
and obtained solutions for the vector ^ defined by Eq. (3.30). The func- 
tion 2 . is sn implicit function of the angular velocities or ^ . In 
the neat part of our work we will discuss a few other features of the system 
before discussing the solutions for ^ or gi . 
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